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TECHNICAL PAPER 


LIMIT CYCLE VIBRATIONS IN TURBOMACHINERY 
I. INTRODUCTION 


High speed turbomachinery is used in many industrial and government applications. 
These include natural gas compressors, turbines in electric power plants, jet engines, and rocket 
engine turbopumps. Despite the differences in application, these machines are all susceptible to a 
variety of common vibration problems. The most common vibration problem encountered with any 
rotating machinery is rotor unbalance vibration. This problem has received much attention over 
the last century and is well understood. Another common problem is dynamic instability. This 
problem is well understood for situations that can be treated using linear analysis techniques; 
however, the limit cycle instabilities that can occur in nonlinear systems have not received as 
much attention. A third problem encountered in high speed turbomachinery is subharmonic 
response to unbalance excitation. This phenomenon requires either nonlinear system character- 
istics or time varying coefficients. This problem has been studied extensively for simple nonlinear 
systems. A significant amount of study has been done for turbomachinery experiencing this 
problem. 

These last two problem areas have each received adequate treatment individually. 
However, situations have occurred where a turbomachine was found to be susceptible to both 
limit cycle instability and subharmonic response for apparently the same operating conditions. 
For example, the high pressure fuel turbopump (HPFTP) of NASA’s space shuttle main engine 
(SSME) has a history of exhibiting subsynchronous vibration at frequencies ranging from 47 to 
56 percent of shaft rotational speed. This includes many occurrences at exactly 50 percent which 
could be attributed to a limit cycle instability or a subharmonic response. This appears to be a 
case of what Hayashi 1 refers to as transition between almost periodic oscillations (limit cycle) 
and entrained subharmonic oscillations. This makes it difficult to determine whether a 
(catastrophic) divergent instability is impending or whether the machine is experiencing sub- 
harmonic response. An examination of this condition is the focus of this report. Before proceeding 
with the discussion of the combination of these two phenomena, a review of the prior work for 
each individual problem is in order. 

General discussions of subharmonic resonance can be found in many vibration texts such 
as the works of Den Hartog 2 and Timoshenko. 3 More detailed analysis can be found in nonlinear 
vibration texts such as that of Hayashi. 1 These works deal primarily with a single nonlinear 
second order equation representing 1 degree-of-freedom (DOF). Tondl 4 provides a thorough 
treatment for 2-DOF models of rotating machinery. Ehrich 5 and Childs 6 have each published 
analyses of rotating machinery having nonlinearity in the form of a clearance or deadband in the 
restoring force. Asymmetry in the nonlinear restoring force was required to demonstrate sub- 
harmonic response. This asymmetry will occur in the presence of deadband when a static load is 
present or when rotor-stator misalignment exists. Bently 7 performed an experimental study 
using a laboratory model rotor. Each of these researchers reached essentially the same conclu- 
sions regarding subharmonic response. Namely, the resonance of a nonlinear rotordynamic sys- 
tem can be excited by a frequency that is near an integer multiple of it. The range of excitation 


frequencies for which this phenomenon will occur depends on the nature of the nonlinearity and 
the other system characteristics. The response will be such that the nonlinear resonance fre- 
quency is tuned to a fraction of the excitation frequency (e.g., 1/2). Hayashi 1 described this as 
subharmonic entrainment. 

Thorough treatments of the causes of dynamic instability can be found in many references. 
Vance 8 devotes an entire chapter to the subject and provides an excellent bibliography. Ehrich 9 
provides a general survey of the fundamental causes of instability and gives some insight into 
means for avoiding them. With regard to the type of instabilities considered here, Ehrich states: 

. . the unifying generality is the generation of a tangential force, normal to an 
arbitrary radial deflection of a rotating shaft, whose magnitude is proportional to (or 
varies monotonically with) that deflection. At some 'onset' rotational speed, such a 
force system will overcome the stabilizing external damping forces which are 
generally present, and induce a whirling motion of ever-increasing amplitude, 
limited only by nonlinearities which ultimately limit deflections." 

The tangential force can be generated by a variety of sources. Among these are fluid bearings and 
seals, turbine aerodynamic forces, and shaft internal damping forces. These tangential forces, as 
well as the restoring and dissipative forces, are usually represented by linear stiffness and 
damping coefficients. The coefficients for the tangential force are usually called cross-coupled 
stiffness and damping, respectively. For fluid bearings and seals, the cross-coupled coefficients 
are directly related to the direct coefficients due to the fundamental physics involved. Black 10 and 
Muszynska 11 both point out that the relationship is due to the transformation from the coordinate 
system which rotates with the fluid average velocity to the inertial system. This fluid average 
velocity is usually slightly less than half the rotor surface speed. This gives rise to the occurrence 
of instability at near one-half rotor speed. A similar relationship is true for internal damping 
except that the coordinate system rotates at the rotor speed. The onset speed of instability is 
determined by examining the stability for various speeds and observing the value of speed that 
causes the rotor to be unstable. This can be accomplished using Routh-Hurwitz techniques or by 
calculating the system eigenvalues. For fluid bearings and seals, it can be shown that the onset 
speed is approximately equal to the system resonance divided by the ratio of average fluid 
velocity to rotor surface speed. Hence, for half synchronous (rotor speed) instability, the onset 
speed is twice the first resonance. 

The restoring, dissipative, and tangential forces discussed above are, in general, nonlinear 
functions of displacement, velocity, and acceleration. Linearizing these functions about an equi- 
librium point as discussed by Gunter 12 results in the linear coefficients used in the stability anal- 
ysis. The onset speed of instability determined by the subsequent analysis addresses stability in 
the small. Due to the nonlinearities, a system may be globally stable when the linearized analy- 
sis predicts an instability. A system in this condition could exhibit a limit cycle instability. 
Sometimes the nonlinearity is in the form of a clearance or deadband in the restoring force. This 
can occur for example when a seal rotor rubs on the stator or when rolling element bearings are 
mounted with clearance between the outer race and the support. The nonlinear restoring force can 
be represented by a piecewise linear function whose value is zero before the clearance is passed 
and whose slope is equal to the "linear" stiffness of the bearing after it is passed. If all other 
force elements are linear, the global onset speed of instability can be determined by assuming the 
value of the clearance to be zero. 13 In this case, limit cycle instabilities can occur at speeds 
below the global onset speed. The frequency of the limit cycle is determined by the nature of the 
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tangential force. For a system with linear fluid forces and clearance in the restoring force, if the 
onset speed of instability is twice the linear resonance, the limit cycle will be at a frequency half 
the rotor synchronous frequency. This is due to the fluid average velocity being approximately 
half the rotor surface speed. If the fluid average tangential velocity is increased or decreased, the 
limit cycle frequency will correspondingly increase or decrease. 

Several researchers have demonstrated limit cycle instabilities using numerical simula- 
tions. Control Dynamics 13 and Day 14 both modeled a single mass rotor with linear direct damp- 
ing, cross-coupled stiffness, and deadband type nonlinear direct stiffness. Control Dynamics 
demonstrated cases that exhibited limit cycles and also cases that did not. They also investi- 
gated the stability in the small about an equilibrium point and concluded that global stability is 
not affected by deadband (clearance). Day also demonstrated cases that exhibited limit cycle 
instabilities and cases that did not. He searched (unsuccessfully) for analytical expressions 
defining the transition points between cases with limit cycle instabilities and cases which only 
exhibit synchronous response to unbalance excitation. Muszynska 15 numerically and experi- 
mentally demonstrated limit cycle instability in which more than one mode (resonance) was 
unstable. As speed was increased, the first mode became unstable. As speed was further 
increased, the second mode became unstable and the limit cycle of the first mode was 
suppressed. 

The understanding of subharmonic response and limit cycle instability as independent 
phenomena is important for many problems that occur in practice. However, of equal importance 
is an understanding of the relationship between each and the response of machinery that is 
simultaneously susceptible to both. High performance turbomachinery falls into this category. In 
fact, any rotating machinery that is operating above a critical speed, has nonlinear restoring 
forces, has static loads and/or misalignments, and has tangential fluid forces can be susceptible 
to both phenomenon. When the frequency of the limit cycle instability is close to a fraction of rotor 
speed (e.g., 1/2) it may be confused with subharmonic response. In fact, if the results of Hayashi 1 
extend to more complex systems, the limit cycle can become entrained by subharmonic response. 
A relationship between the two can be intuitively expected. In simple terms, a subharmonic 
response is caused by the (response amplitude dependent) frequency of a nonlinear system's 
resonance tuning itself to be at a fraction of the excitation frequency so that it will be reinforced. 
Likewise, the limit cycle instability is due to the frequency of a nonlinear system's resonance 
tuning itself so that the (response amplitude dependent) real part of the eigenvalue is zero, 
yielding a sustained transient. 

No other work is known to the author which explores the simultaneous susceptibility to 
these two phenomena, or the relationship between them in rotordynamic systems. The previous 
work on subharmonic response generally neglected self-excitation forces that could lead to 
instability. The previous work on instability did not include the proper conditions to cause sub- 
harmonic response. The work by Control Dynamics 13 was contrived to represent fluid forces for a 
seal with average velocity exactly half the rotor surface velocity. This resulted in a limit cycle at 
exactly half rotor synchronous frequency. It is impossible to differentiate between the two in this 
case. 


Hayashi 1 examines the behavior of a self-oscillatory second-order system (van der Pol's 
equation with periodic forcing term) in the transition between almost periodic oscillations (limit 
cycle whose frequency is not a rational fraction of the excitation frequency) and subharmonic 
entrainment of these limit cycles. He defines (in terms of the external force parameters) regions 
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in which entrainment will occur. A similar determination is proposed as part of this work; how- 
ever, due to the nature of the nonlinearity to be studied here (deadband in the radial restoring 
force) a purely analytical treatment becomes intractable. A combination of numerical and analyti- 
cal approaches must be used. In addition, since it will be desirable to investigate complex real- 
istic systems, a comprehensive numerical analysis tool is needed to model and simulate the 
systems. Such a tool has been developed by the author in order to conduct the proposed 
research. The tool consists of a package of computer programs to perform linear stability, linear 
harmonic response, and nonlinear transient analyses of general turbomachinery. 

The objective of this research is threefold: (1) to characterize limit cycle instability and 
subharmonic entrainment and determine interrelationships between them, (2) to determine 
regions in parameter space for the existence of each and thereby establish criteria for their avoid- 
ance, and (3) to attempt to provide guidance for the interpretation of test data with regard to 
impending divergent instabilities based on observation of subharmonic response or limit cycle 
instability. The investigation begins using the nondimensionalized equations of a single mass 
rotor with the appropriate characteristics. This study employs both analytical and numerical 
techniques. The extension of the single mass model results to a complex, realistic system is 
demonstrated by examining the HPFTP of the SSME. Available test data are examined and 
linear analyses and nonlinear simulations performed. 


n. MODEL FORMULATION 


The initial model used in this study is a greatly simplified representation of a turbopump. 
The model possesses only 2 DOF and, yet, it contains all the characteristics that are germane to 
the phenomena being studied. A model which provides a more complex representation of a 
turbopump is developed in appendix A. The simplified model can be obtained from the more 
complex one, however, its development will be included here for clarity. 

The simplified turbopump model consists of a single mass supported on symmetric 
supports. The supports represent rolling element bearings or fluid film bearings and seals. They 
are treated as linear spring and damper elements with the exception of including clearance or 
deadband in the rolling element bearing force deflection relationship. The shaft flexibility would 
also be included in the linear support spring. Fluid film bearings or seals require the inclusion of 
cross-coupled stiffness and damping terms in order to characterize their influence on rotor 
behavior. The model is excited by three different sources. Rotor mass unbalance provides excita- 
tion at the shaft spin frequency (synchronous). Circumferential pressure distributions in a 
turbopump are represented by fixed direction loads (side loads) at zero frequency. Random noise 
is used to represent a variety of broadband random excitations that may exist in a turbopump and 
to serve as a perturbation to investigate the behavior of the nonlinear system. Figure 1 provides 
a schematic diagram of the system and defines the coordinate systems used in the derivation of 
the equations of motion. 
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Figure 1. Schematic of simplified rotor model. 


The equations of motion are derived from a straightforward application of Newton's 

A A 

second law of motion. Application of this law along the i and j axes yields 


ma y = X Fy , (1) 

ma z = E F z , (2) 

and along the ? r and ? { axes yields 

ma r = ^F r , (3) 

ma, -Yu F t . (4) 

The component forces can be categorized as nonlinear support forces, linear support 
forces, and excitation forces. The force-deflection curve for the nonlinear support force is shown 
in figure 2. The displacement vector diagram is shown in figure 3a. The force in the direction of 
radial displacement R can be written directly as 

F = l-Jt„(R-S) if R>8 

rn I 0 if R<5 • (5 ) 


A A 

This radial force can be resolved into its components along the i and j axes yielding (for R>8 ) 


( 6 ) 
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( 7 ) 



The linear support forces can be written in radial and tangential component form as 


F n er+F tl ?, = (- kiR-cR ) e r +(-cR0+QR) 7, . ( 8 ) 

Transforming these into cartesian components yields 

F yi i + FJ = (- k t Y- cY- QZ ) i+ (- k t Z - cZ+ QY)j . (9) 

The fixed direction side loads and random noise forces are shown in figure 3b. The unbalance 
excitation force is a natural consequence of the proper differentiation (twice) of the displacement 
vector of the center of mass. This vector can be written in either cartesian or polar coordinates as 
the sum of the displacement of the geometric center and the mass eccentricity vector. 

R cg = (R+e cos ( (5-6))7 r +e sin ((5-6)7, 


= (Y+e cos ((5))i + (Z + e sin (f5))j . 

The angle (5 is defined as the rotation angle of the rotor and can be written as 


(10) 


i 3=\ [ a(t*)dt*dt = [ co(t)dt , 

JO JO Jo 


(ID 


where a is the angular acceleration and co is the angular velocity of the rotor. Differentiating 
equation (10) twice with respect to time yields 

( R cg ) = [R-6 2 R-f5 2 e cos ((5 -6) -fie sin ((5-d)\e r 
dt 2 


+ (R6+2Rd-fi 2 e sin (f5-6)+fie cos ((5-6)]?, 


= [Y-fi 2 e cos ((5)-fie sin (/})]/ +[Z-fi 2 e sin ((5)+fie cos ((5)]j . (12) 

The appropriate expressions (cartesian or polar) from equations (5) through (12) can be 
substituted into equations (1) and (2) or equations (3) and (4), respectively, yielding the equa- 
tions of motion for the system. These equations in cartesian form are 

mY+cY+kiY+QZ+kJl- 
mZ+cZ+kiZ-QY+k n il- 


u(R-8)Y = mefi 2 cos ((5)+mefi sin (f5)+F sy +F ny , 
R 


A 

R 


u(R-8)Z = mefi 2 sin ((5)-mefi cos ( p)+F sz +F nz , 


(13) 

(14) 
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and in polar form are 


mR+cR+kiR-m9 2 R+k„(R-8)u(R-8) = mefi 2 cos ((3-9)+me(3 sin ((3-9) 
+(F sy +F ny ) cos (9)+(F sz +F nz ) sin (9) , 

mR9+2mR9+cR9-QR = me ft 2 sin ((3-9)-me(3 cos (f3-9)-(F sy +F ny ) sin (9) 

+ (Fsz+Fnz) cos (9) , 


where 


u(R-8) 


1 if R > 8 
0 if R<8 


(15) 


(16) 


(17) 


Noting from equation (11) that [3 = eo and (3 = a and rearranging the linear and nonlinear 
stiffness terms yields 


mY+cY+k 


1-Y(1-u(R-5))-y&u(R-5) 

R 


Y+QZ = meet) 2 cos (p) 


+mea sin ((3)+F sy +F ny , 


(18) 


mZ+cZ+k 


l-Y(l-u(R-8))-Y&u(R-8) 

R 


Z-QY = meet) 2 sin ((3) 


-me a cos (P)+F sz +F„ z , 


(19) 


in cartesian coordinates and 


mR+cR-m9 R+k 


\-Y(1-u(R-8))~y «(/?-5)j R = meco 2 cos ((3-9)+mea sin ( (3-9 ) 


^(F sy+F ny ) cos (9)+(F sz +F nz ) sin (9) , 

mR9+2mR9+cR9-QR = meet) 2 sin ((3-9)-mea cos ((3-9) 
~ (F S y+F ny ) sin (9)+(F sz + F nz ) cos (9) , 

in polar coordinates where 

k ~ ki "I - kfi y 

and 


y — kjl_ 

1 k 


( 20 ) 


( 21 ) 


( 22 ) 


( 23 ) 
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The bracketed expression multiplying k in equations (18) to (20) is arranged in a way that high- 
lights the two regions of this nonlinear function. For R < S, the third term within the brackets is 
zero and the expression reduces to the linear form [ 1 — y ] . For R > 5, the second term is zero and 


the expression reduces to 


1 -ys. 

. 1 /?J 


In order to generalize the results from this study, it is advantageous to reduce these 
equations to a dimensionless form. This also has the benefit of reducing the number of parame- 
ters in the problem. The first step is to divide equations (18) through (21) by m. When this is 
performed, it is convenient to make the following definitions: 


-L=0) 2 , 

KM 7 


m 


£=2 


m 


(24) 


(25) 


a n and £ correspond to the undamped natural frequency and damping ratio for this system with 
the deadband (<5) and the cross-coupled stiffness (Q) both set to zero. Substituting these defini- 
tions yields the cartesian equations 


Y+2Ca> n Y+Q)£ 


1 -y ( l-u(R-5) )-y & u(R-8) 
R 


+ea sin (/J)+^ + ^ 
^ m m 


Y+ — Z= ew 2 cos (p) 
m ^ 


(26) 


Q 


Z+2£oo n Z+ (o£ l-y(l-u(R-8))-r- Q -u(R-S) Z- ^ Y = ecu 2 sin (p) 
L R \ m 


-ecc cos (P)+^ + ^ , 

^ m m 

and the polar equations 

R+2Co) n R-d 2 R+co^\l-Y(l-u(R-5))-Y^-u(R-5)]R = eco 2 cos (p-6)+eccs in (p-6) 

L R J 

I cos (e)+( F ^+ F "z )sin (0) , 


(27) 


F sy+ F ny 


(28) 


R6+2R6+2Co) n Re-^R = eco 2 sin (p-6) - ecc cos (p-6) 

- sin «» <»> • 

Next, time is normalized by the following substitutions: 


(29) 


(O n ’ 


( 30 ) 


9 


( 31 ) 


r d% r. 

(32) 

In these expressions, | represents any variable which is being differentiated. Performing these 
substitutions and dividing through each equation by co% yields the cartesian equations 


r"+2£r+ 1-y ( \-u(R-8))-y£- H(fl-<5) Y+qZ= ep 2 cos (p) + ep sin (p) + -^- + 


F sy F n y 


k k (33) 


Z"+2£Z'+[l-y ( l-u(R-S ) )-/-& u(R-8)] Z-qY = ep 2 sin (p) - ep cos (/J) + ^ ^ , 

R J k k (34) 


and the polar equations 


R” +2£R' -d ,2 R + \\-y {\-u{R-8 ))-y u(R-8)]r = ep 2 cos (p-0)+ep sin ( p-6 ) 

L R J 


F +F 

1 sy ^ 1 ny 


cos (6) + sin ( 0 ) , 


R&'+2R'ff+2ZRff-qR = ep 2 sm (p-d)-ep cos (ft-fl) - T Jy+Fny sin ( 6 ) 

' k ' 


F„+F l 


^■J cos (0) , 


where 


= _Q_ = Q 

^ m(D% k 


P = - GL , 

F c On 


(On 


The angle P defined by equation (11) can be expressed in terms of x, p, and p as follows: 


P= f f 6)%P(T*) = I f p(X*)dX*dX = [ 1 p(x)dx 
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The final step in the nondimensionalization process is to define the following new vari- 
ables 

Y=y8, (41) 

Z = z8, (42) 

R = r8 . (43) 

Substituting equations (41) through (43) into equations (33) through (36) and dividing through 

by 8 yields the cartesian equations 

y" +2£y' +[l-y (l-w(r-l)) u(r-l)]y+qz = ap 2 cos (p) + ap sin (/ 3)+g y +ri y , (44) 

z"+2£z'+[l-y(l-u(r-\)) u(r-\)\z-qy = ap 2 sin (p)-ap cos (p) +g z +ri z , (45) 

and the polar equations 

/'+2£r'-0' 2 r+[l-y (l-M(r-l)) -y u(r-l)] r = ap 2 cos (p-6)+ap sin (p-6) 

+ ( gy+riy ) cos (0) + (g z +Vz) sin (6) , (46) 

re"+2r'0'+2Crd'-qr = ap 2 sin (p-6) - ap cos ( p-6 ) - ( g y +ri y ) sin (0) 


where 


+ (g z +r lz) cos (0) . 



gz = 

n y = 

Vz = 


r sy 

Ts ’ 

Fs± 
k8 ’ 

F* ny 

Ta ’ 

k8 


(47) 

(48) 

(49) 

(50) 

(51) 

(52) 


These equations contain a total of 10 dimensionless parameters. The number of parameters can 
be reduced by making a few assumptions. First, since the system is symmetric, g z can be 
assumed to be zero without loss of generality. Second, for the cases to be investigated here, 
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p will be very small and can be neglected. Finally, rj y and Tfc are assumed to be uniform random 
number sequences with the same range These assumptions yield a total of seven 

dimensionless parameters: g y , a, p, q, y, and rf. A complete parametric study including all 
seven parameters is not feasible; however, a subset of the most significant will be selected 
based on the analyses discussed in the next section. 


III. ANALYTICAL TREATMENT 


The system model defined in the previous section can be analyzed in several ways. The 
equations can be linearized about various equilibrium conditions in either the cartesian or polar 
form. This approach is appropriate when examining limit cycle instability. The forced response 
can be determined for the nonlinear system using a harmonic balance procedure. This is appro- 
priate when examining subharmonic resonance. This method requires a combined analyti- 
cal/numerical approach in order to determine a solution. Each of the linearizations and the har- 
monic balance method yield certain insights into the characteristics of the system; however, no 
single approach is adequate to describe the general case. 


A. Linear Stability 

The simplest approach to linearization is to neglect the deadband (8) by assuming it to be 
zero. For this approach, equations (33) through (36) must be used since equations (44) through 
(47) are based on nondimensionalizing equations (33) through (36) using 8 as the basis. As 
discussed in reference 13, the stability determined in this manner is the global stability. Working 
with the cartesian equations, the Routh criterion will be applied to the characteristic equation for 
the system to determine the stability boundary. Applying the 8 = 0 assumption to equations (33) 
through (34) yields 


Y"+2CY'+Y+qZ= ep 2 cos (p) + epsin (p) + ^ + ^ , 

K K \J J ) 

Z"+2CZ'+Z-qY= ep 2 sin (/?) - ep cos (p) + Ess. + Esl ^ 

k k (54) 

Stability of the system is determined from the homogeneous solution to these equations. 
Assuming solutions of the form 


Y=Ye Xl , 
Z = Ze Xt , 


results in the following algebraic equation: 

A 2 +2£A+1 q 
-q A 2 +2£A+1 



(55) 

(56) 


(57) 
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This equation has a nontrivial solution only when the determinant of the matrix is zero. 
Expanding this determinant yields the characteristic equation for the system, 

A 4 +4£A 3 + (4£ 2 +2)A 2 +4£A+l+<? 2 = 0 . (58) 


Applying the Routh test to this polynomial yields the following array: 


A 4 

1 

A 3 

4C 

A 2 

4C 2 + 1 

A 1 

4C(4C 2 -? 2 ) 

4£ 2 + l 

A« 

l+q 2 


4^+2 l+q 2 

4C 

l+q 2 


For positive £, the only possible sign change in the first column occurs in the A 1 row. This term 
will remain positive as long as 



(59) 


Therefore, equation (59) is the stability criterion for the system. In terms of the original system 
parameters given in equations (18) through (21), this condition becomes 



(60) 


As discussed in section I, there is generally a relationship between Q and c (or q and 2£) that is 
determined by the fluid dynamics of the particular system under consideration. This relationship 
is generally a function of rotor speed. The form of this relationship that is assumed in this study 
is discussed later. 

It is of interest to determine the roots of the characteristic equation (equation (58)) in the 
marginally stable condition (-$-= l] . Equation (58) can be expressed in factored form as 


(A 2 +2£A+1 +jq) (X 2 +2£X+l-jq) = 0 . (61) 

The roots of this equation can be determined by applying the quadratic formula to each term 
yielding 


A,- = -f±Vf 2 -l ±jq . 

Substituting q = 2£ into equation (62) yields the four roots 

Ai,2 = ±jl , 


(62) 


(63) 
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A 3 ,4 = -2C±/1 . 


(64) 


These roots can be related back to unnormalized time yielding 

*1.2 = ± j ( °n , (65) 

A 3 ,4 = -2 CcOn±jco„ . (66) 

It can be seen from comparing equations (59) through (66) that the frequency of the instability is 
equal to (o n and the instability occurs when the ratio y is equal to 0)„ . It will be shown next that 

for the nonlinear system (5* 0) the frequency of the limit cycle that can occur when y is less 
than a )„ is equal to ^ for certain cases. 

B. Self-Excited Equilibrium 

The approach taken next is to seek an equilibrium solution to the nonlinear, homogeneous 
equations in polar form. The homogeneous form of equations (46) and (47) is first rewritten for 
two cases, r < 1 and r > 1. For r < 1 the result is 


r'+2£r-d' 2 r+(l-Y)r = 0 , 

(67) 

rd"+2r'$'+2£r6'-qr = 0 . 

(68) 


For r > 1 the e t equation is identical to equation (68) and the e r equation is given by 


r"+2£r'-d' 2 r+r-y= 0 . 


(69) 


For an equilibrium solution, the quantities r", r\ and 6" are all zero. Eliminating these terms from 
equation (68) yields the equation governing the equilibrium angular velocity 


2£r o d'o-qro = 0 , 


or expressed explicitly 



(70) 


(71) 


This applies for all values of r. It should be noted that this is the ratio that governs the system 
global stability as shown in the preceding paragraphs. For the case where r < 1 , eliminating the 
same quantities as before in equation (67) yields the equation governing the equilibrium radius 


(1 — Y ~ & o) r o — 0 . (72) 

The only equilibrium solution to this equation is r = 0 which was expected since, for this range of 
r, the system is linear. For the case where r > 1, the same quantities are eliminated from equa- 
tion (69) yielding 
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(73) 


d-di)r 0 -Y=0 , 


or explicitly 

r 0 = - 
1 


In order for a nonzero solution to exist for ro, ro must have a value greater than 1. Examination of 

q 

equation (74) shows that to meet this condition, ^ must be greater than Yl-y . This can be 

recognized as the instability criterion for the linear system obtained for the case r < 1. In other 
words, if the linear system defined for the case of the radial deflection being less than the dead- 
band is stable, no equilibrium solution (other than ro = 0) is possible. If the linear system is 
unstable, the equilibrium radius is governed by equation (74). This equilibrium motion represents 
a limit cycle instability. It can also be seen from equation (74) that at the global stability 
q_ 

boundary, 2 £= 1, the solution is unbounded, as would be expected. 

C. Mass Unbalance Equilibrium 

The self-excited equilibrium determined in the above paragraph gives insight into the 
behavior of the system; however, it does not represent a very practical case since it is completely 
unforced. Another case of interest is the equilibrium solution to the system when excited only by 
mass unbalance. The solution is assumed to be harmonic with frequency equal to the excitation 
frequency (synchronous response), i.e., 6' o = p.The equilibrium assumptions are applied to 
equations (46) and (47) (assuming p = 0). The result for r < 1 is 

(1-y-p 2 ) r 0 = ap 2 cos (fi-6) , ( 75 ) 

(2 £p-q) r 0 = ap 2 sin ( p-6 ) . (76) 


7 _ 

a '2 
0 


1- 


2C 


(74) 


For r > 1 the 7, equation is identical to equation (76) and the 7 r equation is 

(1-p 2 ) r 0 -y= ap 2 cos (P-6) . 

Equations (75) and (76) can be solved simultaneously for ro yielding 

,2 


r o = - 


a P 


V( l-y-p 2 ) 2 + (2£p-tf) 2 

for r < 1. Equations (77) and (76) can be solved simultaneously for r 0 yielding 


(l-p 2 )7±Va^p 4 (l-p 2 ) 2 -i-a 2 p 4 (2^p^^) 2 -y 2 (2^p-^) 2 

(l-p 2 ) 2 +(2Cp-<?) 2 


(77) 


(78) 


(79) 
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for r > 1. Equation (78) is valid for all values of the parameters for which the linear system 
defined for the case of r < 1 is stable. Three conditions must be met for equation (79) to be valid. 
First, the quantity under the radical must be greater than zero since r must be real valued. 
Second, one of the differential equations from which equation (79) was derived (equation (69)) is 
only valid for r > 1; therefore, the predicted value of r must be greater than 1. Third, the equilib- 
rium must be stable (determined in a following section). 


D. Side-Force Equilibrium 


A third equilibrium point of interest is the case with fixed direction side forces and no 
unbalance. The solution is assumed to be static, i.e., all time derivatives are equal to zero. 
Applying these assumptions to equations (46) and (47) (and recalling that g z was assumed to be 
zero in section II) yields 

(l-y)^o = gy cos (6q) , 


qr 0 = g y sin (0 O ) , 


(81) 


for r < 1. For r > 1, the ?, equation is identical to equation (81) and the e r equation is 

ro~Y= 8y cos (6 0 ) . 

Solving equations (80) and (81) for the case of r < 1 yields 

8y 

V(l-y ) 2 +$ 2 
60 = lan " 1 (ir^) • 

This solution can easily be expressed in terms of cartesian coordinates as 


(82) 


(83) 

(84) 


( i - y ) 2 +<? 2 

. _ m 

zo = . 

(l-y) 2 +4 2 

For the case of r > 1, equations (81) and (82) yield the following quadratic equation in rg. 


(85) 

( 86 ) 


(l+4 2 )ro-2yr o +y 2 -g 2 = 0 . 
The solution of this equation is found to be 


(87) 
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( 88 ) 


r o = 


7+V gyd+g 2 )-g 2 y 2 

iV 


It can be shown that for all cases where ro predicted by equation (83) is greater than 1, the 
quantity under the radical in equation (88) is greater than zero. The corresponding solution for 0o 
is 


do = tan 1 


I m \ 

U-yJ 


(89) 


Equations (88) and (89) can be transformed into cartesian coordinates yielding 

y 0 = r 0 cos (d 0 ) , (90) 

zo = r 0 sin (0 O ) • (91) 


E. Linearization About Self-Excited Equilibrium 

Equilibrium solutions have now been determined for three specific cases of interest. 
These cases are the self-excited system (equations (71) and (74)), the synchronous response 
to unbalance excitation (equations (78) and (79)), and the static response to side forces 
(equations (83) through (91)). Linearizations of the equations of motion can be determined for 
each of these cases. The polar form of the equations will be used for the first two cases and the 
cartesian form for the last. 

Linearization of the equations of motion is only necessary for r > 1 since the system is 
already linear for r < 1. The linearized form of the equations of motion for the first two cases can 
be developed together. The appropriate excitation terms can be dropped in order to obtain the 
equations for the self-excited case. Equations (46) and (47) are used with only the unbalance 
excitation forces retained. Also, as in the determination of the equilibrium solution, the angular 
acceleration fi is assumed to be zero. The resulting equations are 

r"+2£V-0' 2 r+r-y= ap 2 cos (/J-0) , (92) 

rd"+2rd'+2£rd'-qr = ap 2 sin (/J-0) . (93) 

The linear equations will be obtained by examining perturbations about the equilibrium solutions. 
This is achieved by making the following substitutions into equations (92) and (93): 


r = ro+r , 

(94) 

d' = d' 0 +d ' . 

(95) 
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The variables r and 6' are the perturbation variables. The linearization is accomplished by 
making the substitutions and neglecting terms of second order or greater in the perturbation 
variables. 


For the self-excited case the resulting equations are 

r"+2£r'+(l-6'o)r-26 , oro9 , = r 0 (6'o-l)+Y , (96) 

rod"+2d'or'+(2£d'o-q)r+2£rod' = (« q-2CO' 0 )r 0 . ( 97 ) 

The equilibrium conditions are obtained by setting the perturbation terms to zero. The results are 
identical to those obtained previously (equations (71) and (74)). The stability of this equilibrium 
is determined from the homogeneous solution of equations (96) and (97). The characteristic 

equation for this system (with the substitution d'o = ql2£) can be written as 


\ 4 C 2 I \ 2 C ) (98) 

It can be shown that the Routh test applied to equation (98) yields the same stability criterion as 
the global stability criterion given in equation (59). 


F. Linearization About Mass Unbalance Equilibrium 

For the synchronous response case, the linearized equations are (with the substitution 
0'o = P) 


r"+2£r'+(l-p 2 )r-2prod' = ro(p 2 -l)+y+ap 2 cos ((5-0) , 


(99) 


rod" + 2pr' + (2£,p-q)r+2C,rQ6' = ( q-2£p)ro+ap 2 sin (/ 3-0 ) 


( 100 ) 


Before proceeding with the perturbation analysis, the sine and cosine terms must be expressed 
in terms of the perturbation variables. For the constant speed condition under consideration, it 
can be seen from equation (40) that (5 = pt+fi o- 9 can be written in terms of the perturbation 
variable as 


-f 

Jo 


0= \ (O'o+O 9 )dx- pr+0o+ I 0'dr=pT+$o+d 


(V. 

Jo 


( 101 ) 


The resulting argument of the sine and cosine terms in equations (99) and (100) is Pq-0o-0. 
The equilibrium solution for ro obtained by setting the perturbation terms to zero is identical to 
that given in equation (79). The equilibrium solution for the phase angle fo- &o can be expressed 
implicitly by 
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(102) 


cos ( po-do ) = 


(1-pVo-y 



sin (po-do) = 


(2 £p-q)r 0 



(103) 


The terms cos (po-do-d) and sin (p 0 -d 0 -d) can be expanded using elementary trigonometric 
addition identities into the form 


cos ( po-do-d ) = cos (/)o-0o) + sin (Po~do)d , (104) 

sin (po-do-d) = sin ( po-d 0 )-cos ( Po-d 0 )d . (105) 

The small angle assumptions cos (d) = 1 and sin (0) = d have already been incorporated into 
these equations. Equations (102) and (103) are now substituted into equations (104) and (105). 
The resulting expressions are substituted into equations (99) and (100). Retaining only those 
terms involving the perturbation variables yields 


r"+2£r'+(l-p 2 )r-2prod'-(2£p-q)rod =0 


(106) 


r 0 d"+2p7+(2Cp-q)?+2Crod'+((l-p 2 )ro-Y)d = 0 . 

The characteristic equation which results from these homogeneous equations is 
A 4 +4CA 3 +(4C 2 + 2p2+2-^)A 2 + 2[2C(l-p 2 )-^+2p(2Cp-9)]A 

+(2Cp-?) 2 +(i-p 2 )(i-p 2 -^)=o . 


(107) 


(108) 


This equation involves the equilibrium solution ro and, hence, a stability criterion cannot be easily 
expressed analytically. The stability of this equilibrium will be determined numerically for various 
values of the system parameters and presented in section IV. 


G. Linearization About Side-Force Equilibrium 


The side-force equilibrium case will now be treated using the cartesian form of the equa- 
tions of motion. The equations of motion (equations (44) and (45)) with only side-force excita- 
tion can be written as 


y"+2£y'+qz-g y 


y 

Vy 2 +z 


y . 


z"+2£z'-qy-g z 



(109) 

( 110 ) 
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where the nonlinear restoring force has been isolated on the right-hand side and r has been 
expressed explicitly in terms of y and z. The right-hand sides of equations (109) and (110) can 
be linearized using their Taylor's series expansions about the equilibrium points yo» zo- Retaining 
only the first-order terms and rearranging the equations results in the homogeneous equations 


r+2cr+ 1 - 


+ - 


ypy 


yozoY 

Vyo+zo (yo+zo)*l l(>o +z$)*l 


y + i- 


\ z + qz = 0 


(111) 


z"+2£z'+l 1 


V 


■ + 


yo+z6 


i |r | 

(yo+zo)*} 


z + 


r yozoY ' 
-(yo + *o)'J 


y-qy = o . 


(112) 


The equilibrium points (given in equations (90) and (91)) were derived for the case where g z = 0. 
In examining equations (111) and (112), it would be convenient to be able to consider an equilib- 
rium point on one of the coordinate axes, such as zo = 0. This can be accomplished without loss of 
generality by considering the load to be applied at an angle of - do with respect to the y axis 
where do is given by equation (89). In this condition, yo will be equal to ro as defined by equation 
(88) and zo will be zero. Since the system is symmetric, the stability of this equilibrium point will 
be identical to the original point. For this case, equations (111) and (112) become 

y”+2£y'+y + qz = Q , 
z" + 2£z'+jl -j^jz-qy = 0 . 


(113) 

(114) 


The corresponding characteristic equation is 

A 4 + 4fA 3 + (2-^ + 4f 2 )A 2 +2f(2-X)A + l-^ + ,2 = 0 . 


( 115 ) 


Application of the Routh test to this equation yields the following stability criterion: 

(116) 

Since this relation involves the equilibrium solution tq. A general conclusion about the stability 
cannot be drawn. However, for cases where the quantity within the parentheses is positive, this 
requirement is less restrictive than the global stability criterion (equation (59)). This implies that 
over some range, the addition of side forces to a system with deadband in the restoring force has 
a stabilizing influence (in the small). More general results will be presented in section IV by 
evaluating the roots for various values of the system parameters. 


H. Subharmonic Response Analysis 

The stability of the system model has been examined for various equilibrium conditions. 
The solutions obtained provide insight into the characteristics of the limit cycle instability. An 
understanding of the characteristics of the subharmonic response will now be sought using a 
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harmonic balance method. The dimensionless cartesian equations (equations (44) and (45)) will 
be used. The assumptions made previously (g z = 0 and fi = 0) will be made here, and the random 
noise excitation (rf y and 7] z ) will be neglected. 

The nonlinear restoring forces in equations (44) and (45) can be rewritten to facilitate the 
development of the harmonic balance equations. The modified equations are 

y"+2£y'+(l-y)y+qz+f y (y,z) = ap 2 cos (pz)+g y , (in) 

z"+2Cz'+(l-Y)z-qy+f z (y,z ) = ap 2 sin (pz) , ( 118 ) 


where 


and 


f y (y>z) 


y 1 1 — = — \ y if Vy 2 +z 2 > 1 

I Vy 2 +z 2 / 

0 if Vy 2 +z 2 < 1 



1 \ z if Vy 2 +z 2 > 1 

Vy 2 +z 2 / 

0 if Vy 2 +z 2 ^ 1 


(119) 


( 120 ) 


The method used here is essentially the same as that used by Noah. 16 This investigation will be 
limited to subharmonics of order 2 (one-half synchronous). The solutions for y and z are assumed 
to be superpositions of a fundamental sinusoidal component and N of its harmonics. The funda- 
mental in this case is the one-half synchronous subharmonic. These solutions are given by 


N 

1 

n= 1 

N 

x 

n= 1 


Ur)- 

- K. sin 



1 2 I 


\ 2 >1 

(121) 

U rl 

- sin 

K r ))- 


v 2 } 


\ 2 II 

(122) 


The nonlinear restoring forces f y and f z are approximated by similar harmonic expansions given 
by 


N 

fy = c yo + cos 

n= 1 

( n 2 r )~ dy * sin 1 

K r ))- 

(123) 

N ( 




fz = Qo + Z [ C Z» cos 

( n o' T ) ” sin 1 

I n 2 T )) • 


«= 1 

v Z f 


(124) 
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Substituting equations (121) through (124) into equations (117) and (118) and performing a 
harmonic balance yields a set of AN +2 linear equations in 8W+4 unknowns ( a yo , a Zo , c yo , c Zo , a y ., 
a z ., b y ., b Zj , c yi , c Zi , d yr d Zi , .... i = \,N). The additional AN +2 equations needed for a solution are 
determined from the relationship between f y and f z and the solutions y and z given by equations 
(119) and (120). These equations can be solved using an iterative numerical procedure. The 
results will provide insight into the effects of the various system parameters on the subharmonic 
response. 

Analytical solutions have been developed for the model formulated in section II. These 
solutions are expressed in terms of the model parameters. In section IV, numerical values (or 
ranges of values) will be specified to define the model. The numerical values of the analytical 
solutions will then be presented for various values of the model variables. 

IV. MODEL DEFINITION AND ANALYSIS RESULTS 

A simplified single mass model of a turbopump rotor has been formulated in section II. 
Analytical expressions for various equilibria and linearized stability conditions were obtained in 
section III. In this section, numerical values and ranges of values will be specified in order to 
define the model. Using these values, numerical results will be presented for the solutions 
developed in section III. 


A. Model Definition 

The model consists of seven dimensionless parameters ( g y , a, p, £, q, y, and rf). These 
parameters are defined at the end of section II. The values to be chosen for these parameters will 
be based in part on the author's experience and in part on the analytical expressions obtained in 
section III. Initially, the random noise parameter rj will be neglected. This assumption was made 
in the developments of section III. Parametric studies of the effects of Tj will be performed later 
using simulations. 

The nonlinear stiffness to total stiffness ratio (y) is restricted by its definition to range 
from zero to one. Typically, in rocket engine turbomachinery, rolling element bearings provide a 
significant, if not the majority, of the rotor support stiffness. These bearings frequently are 
mounted with clearance between the outer race and the bearing support. This clearance provides 
the deadband 8 discussed in section II. A value for y of 0.75 has been selected to represent a 
typical rotor support situation where clearance mounted bearings provide a majority of the rotor 
support stiffness. 

The dimensionless side load g y and dimensionless unbalance a were normalized by the 
deadband 8. Hence, they should have values on the order of unity to represent cases where the 
rotor is operating in a highly nonlinear fashion. Values much greater than unity will tend to 
obscure the deadband. Values much less will cause operation in the linear range of the function 
defining the rotor support (equation (5)). Nominal values of 1.0 and 0.5 have been assigned to g y 
and a , respectively. 

The shaft angular velocity p is a primary parameter in any investigation of rotating 
machinery, and a wide variety of values will be examined. However, the upper limit of interest is 
the maximum of the global onset speed, the unbalance stability threshold, and the side-force 
stability threshold. 
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As discussed in sections I and III, the dimensionless cross-coupled stiffness q is 
generally related to the fluid damping 2£. This relationship is due to the fluid average tangential 
velocity being a function of rotor surface velocity. If a fixed ratio is assumed, q can be expressed 
in terms of this ratio, the rotor speed, and the damping. This ratio is typically slightly less than 
one half. If the configuration is such that the fluid is entering the fluid seal or bearing with a 
tangential velocity greater than the rotor surface speed, the ratio can be greater than one half. 
Designating this ratio as a, the cross-coupled stiffness can be expressed as 

q = 2&o . (125) 

Since a primary objective of this study is to examine systems which are simultaneously 
susceptible to limit cycle instability (governed by o) and subharmonic resonance, the value of a 
has been selected to be 0.48 (close to but slightly less than one half). This ratio has frequently 
been observed in SSME test data. Values greater than one half will also be examined. The 
damping ratio £ is typically low in rocket engine turbomachinery. A nominal value of 0.10 has 
been assigned to this study. This value is representative of the damping in the HPFTP of the 
SSME. 


Nominal values have been assigned for all parameters of the model. These values are 
summarized in table 1. The results of the analyses of section III will be presented for these 
nominal values. 


Table 1 . Nominal model parameters. 


Parameter 

Value 

9y 

1.0 

a 

0.5 

p 

0.-4.0 

c 

0.1 

<T 

0.48 

7 

0.75 


B. Analysis Results 

The first result developed in section III was the stability criterion for the case of zero 
deadband (equation (59)). Substituting equation (125) into equation (59) results in the criterion 


<i = ^- = 2.08 . 
° 0 . 


(126) 


As discussed in section III, this is the global onset speed of instability for the system. It will be 
designated by p gl . The roots of the characteristic equation (equation (58)) are shown in figure 4 
for various values of £ and <x as p is varied from 0.0 to 4.0. The plot for variations in £ shows the 
roots in the complex plane. As <r is varied, the roots follow identical loci except that the speed 
correspondence is different. For this reason, the critical damping ratio is displayed for the roots 
as cris varied. As expected, the imaginary axis crossing is determined by equation (126). 
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(a) Effects of £ on roots in complex plane. 
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(b) Effects of (Ton critical damping ratio of roots. 

Figure 4. Roots of characteristic equation (equation (58)) as p is varied from 0.0 to 4.0. 
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The next result developed was for the equilibrium solution to the homogeneous polar 
equations. The results showed the possibility of a limit cycle instability with a frequency given by 
equation (71) and radius given by equation (74). Substituting equation (125) into these 
equations and applying nominal values for the parameters yields 


e' 0 = op- 0.48p , 


and 


0.75 

(1-0.48 2 p 2 ) 


(127) 


(128) 


As discussed in section III, this equation is valid only when ro is greater than 1. This will be true 
whenever q! 2£ is greater than V 1-y . Expressing this in terms of the parameter values yields 


„>!E£ = iMI = i.o4 . 

y <* 0.48 


(129) 


Since this is the initial value for which a limit cycle is possible, it will be designated by p ic . 
Equation (128) is plotted in figures 5 and 6 for a range of p from zero to the global stability limit. 
Figure 5 presents the solution for various values of y and figure 6 for various values of o. The 
primary effect of each parameter is to change the range of values of p over which the 
homogeneous equilibrium is possible, y does not alter the upper limit (p g i) but it has a strong 
effect on the lower limit (p/ c ). For y - 0 (the linear problem), no homogeneous equilibrium is 
possible. For y = 1, the equilibrium is possible for all values of p from zero to the global onset 
speed of instability, o affects the lower and the upper limit since both are inversely proportional 
to o. As would be expected from examination of equation (128), the amplitudes increase 
dramatically as the value of p approaches the upper limit. 

The equilibrium solutions for the unbalance mass excitation case are given by equations 
(78) and (79) for values of r 0 < 1 and r 0 > 1, respectively. The equations are subject to the 
validity conditions discussed in section III. Stability of the equilibrium is governed by the roots of 
the characteristic equation (equation (108)). The solutions are plotted in figures 7 through 10. 
The absence of a portion of a curve indicates that the solution is not valid or is unstable in that 
region. The stability condition for the linear case of ro < 1 (p/ c ) and the global stability condition 
(Pgi) are indicated on the figures. The stability threshold for the equilibrium solution is also 
indicated. This value was determined by examining the real parts of the roots of equation (108) 
as the parameters were varied. This threshold is designated by p ue . 

Figures 7 through 10 present the solutions for various values of y, a, £, and o, 
respectively. These curves contain regions where dual solutions are possible for the same value 
of all parameters. One solution is for ro < 1 and the other for ro > 1. There are also regions where 
no mass unbalance equilibrium solutions are possible. These are the regions where the 
equilibrium would reside in the linear range ro < 1, but the system is unstable in this range 
(p > pic), or the solution would reside in the range ro > 1 and the solution is unstable in that 
range (p > p ue ). In these ranges, some combination of the homogeneous limit cycle solution given 
in figures 5 and 6 and the unbalance might be expected, yhas a strong influence on the limits of 
these ranges. As y is increased, both p/ c and p ue (when ro > 1) decrease. It also causes an 
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Figure 9. Mass unbalance equilibrium solution versus p. Effects of £ on solution amplitude. 



Figure 10. Mass unbalance equilibrium solution versus p. Effects of cron 
solution amplitude, pi c , and p„,. 





effective lowering of the system resonance frequency. This would be expected since y softens the 
support in a nonlinear way. The unbalance parameter a does not affect the limits p/ c or p gi . 
However, for the larger values, the solution never drops below 1.0 and, therefore, it exists for all 
values of p up the equilibrium stability threshold p ue • C has little effect on the speed ranges 
where unbalance equilibrium is impossible. Its primary effect is on the amplitude of the solution 
at resonance (as with linear systems), cr directly affects both p/ c and the global stability limit p g j 
(each is inversely proportional to cr). It affects the equilibrium solution stability threshold p ue in a 
manner similar to its effect on p gl . 

The equilibrium solutions for the fixed-direction side forces are given by equations (83) 
through (86) for r 0 < 1.0 and equations (88) through (91) for r 0 > 1.0. These equations are in 
terms of q which should be replaced with its definition given by equation (125). The stability of 
the solution for ro > 1.0 is determined by equation (116) or direct evaluation of the roots of 
equation (115). Again, the definition of q given by equation (125) should be substituted into 
these equations. The roots of equation (1 15) are plotted in figure 11 for various values of £ and cr 
and in figure 12 for various values of yand g y . In both figures, p is varied from 0.0 to 4.0. As in 
figure 4, the roots for variations in a and g y are presented as critical damping ratios. The results 
for all cases have certain similar characteristics. For p = 0.0, the y and z axis equations are 
uncoupled due to the special form assumed for q (equation (125)). The roots occur in complex 
conjugate pairs with all real parts equal to £. One pair of roots has an imaginary part equal to 

V 1.0— C 2 (corresponding to equation (113)) and the other has an imaginary part equal to 

V 1.0- y/r 0 - C 2 (corresponding to equation (114)). As p increases from zero, the equations 
become coupled and the roots approach each other along the line defined by the real part equal to 
£. When the roots meet, they branch symmetrically away from the vertical line given by the real 
part equal to £. One branch becomes more stable while the other moves toward the right half 
plane. 


The parameter £ affects the roots primarily in two ways. First, it defines the real part of 
the roots for values of p prior to the intersection of the roots. Second, since q is a linear function 
of £, increasing it causes the root intersection and imaginary axis crossing to occur for lower 
values of p. Increasing a has a similar effect without shifting the value of the real parts for small 
p. For some small values of a, the roots never intersect for the range of p examined here, y 
affects the imaginary part of the lower root for p = 0. For smaller y, this root moves closer to the 
higher root. For y- 0 (the linear system), the roots are identical. Since the roots are closer 
together their intersection and the imaginary axis crossing occurs for smaller values of p. The 
effect of g y is similar to that of y, only inverted. This would be expected from examination of 
y Y 

equation (115). yand tq always appear as in this equation and, from equation (88), it can be 

seen that ro is almost a linear function of g y . For large values of g y or small values of y equation 
(116) indicates that the side force equilibrium is less stable than the global stability of the linear 
system. However, for increasing values of g y , the stability condition approaches the global 
stability condition. This result was previously shown by Control Dynamics. 12 

The subharmonic resonance solution was developed in section III using a harmonic 
balance procedure. This procedure was implemented using the Newton-Raphson method as 
presented by Noah. 16 The primary interest here is in the values of the subharmonic components 

of the series solution (a y i, b y i, a z i, b z i). The magnitudes of the y axis components (Va 2 , + b yi ) 
are plotted in figures 13 through 17 for various values of g y , a, y, £, and a, respectively. The z 
axis components behave similarly. The absence of a portion of a curve in these figures indicates a 
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(a) Effects of £ on roots in complex plane. 



(b) Effects of c on critical damping ratio of roots. 


Figure 11. Roots of side-force equilibrium characteristic equation (equation (115)) 

as p is varied from 0.0 to 4.0. 
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Figure 17. Subharmonic response solution versus p. Effects of cr on 
magnitude of y-axis component. 

failure to converge to a solution. This failure most likely indicates that the assumed form of the 
solution does not represent a solution that can actually exist for the given set of parameters. On 
the other hand, the achievement of convergence for another set of parameters does not guarantee 
that the solution will take the prescribed form, only that it is possible. This is important to note 
when examining results from simulations which exhibit nonfractional subsynchronous response 
(e.g., limit cycle instability), since these motions cannot be represented by the subharmonic 
solution form. 

The effect of g y on the subharmonic response can be seen from figure 13. As g y is 
increased, the range of occurrence of the subharmonic narrows. The upper limit moves upward 
slightly. The amplitude of the response is not significantly affected at a given value of p as long 
as the value is within the range of occurrence. Within the range of occurrence of the subharmonic, 
the response increases with increasing p. 

The unbalance parameter a affects the subharmonic response in much the same way as g y . 
This can be seen in figure 14. One difference to note is that the range of occurrence initially 
broadens and subsequently narrows as a is increased. Also, the amplitude at a given value of p 
varies a little more with a than with g y . 

The most significant influence comes from the parameter y This is not surprising since 
this parameter gives a measure of the degree of nonlinearity in the system. The effects of y can 
be seen in figure 15. For y= 0, no subharmonic is developed since the system is linear. As y 
increases to 1 (its maximum), the amplitude and the range of occurrence of the subharmonic 
response increase. The upper limit of the range drops somewhat. 
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The damping parameter £ has much less influence than the parameters discussed 
previously. The primary effect is to increase the upper limit of the range of occurrence of the 
subharmonic. This effect is shown in figure 16. Since £ has little effect on the amplitudes, the 
curves for each successive value of £have been shifted vertically for clarity. The cross-coupled 
stiffness parameter <7 also primarily affects the upper limit. Figure 17 shows its effect on the 
subharmonic response. These curves are shifted as in figure 16. For very small values, no 
subharmonic response occurs. After a certain value is exceeded, the range begins to broaden. 
Both the lower limit and the upper limit are affected. The largest change comes as <r approaches 
a value of 0.5. As this value is approached, the upper limit of the range of occurrence and the 
amplitude of the response at the upper limit increase dramatically. This is due to the coalescence 
of the limit cycle instability frequency with the subharmonic response frequency. The upper limit 
and the amplitude decrease as o increases beyond 0.5. These effects are indications of the 
inherent relationship between limit cycle instability and subharmonic response that was 
postulated to exist in section I. 


C. Preliminary Observations 

The results developed in section III and presented in this section apply for certain 
restricted, sometimes nonrealistic conditions. For example, the homogeneous equilibrium 
solution is of little value for predicting actual system response since any real system will 
possess some amount of excitation. Similar statements could be made about the other solutions 
developed. However, each yields some insight into the characteristics of the behavior of the 
system and some general conclusions can be drawn. The conclusions drawn for the restricted 
cases can be extended to apply for cases where one parameter is dominant. For example, if the 
unbalance is much larger than the side force, the characteristics observed for the unbalance 
excitation case would be expected to hold. When this is not the case, other conclusions can be 
drawn based on the results of the restricted cases. 

The homogeneous solution to the polar equations yielded a range of occurrence, a 
frequency, and an amplitude for the limit cycle instability. The analysis of the unbalance response 
equilibrium determined the possibility of a solution in regions where the limit cycle equilibrium is 
also possible. The stability analysis of the side-force equilibrium demonstrated that this 
excitation increased the stability of the system (in the small). It is reasonable to expect that both 
of these effects would tend to reduce the range of occurrence of the limit cycle. The limit cycle 
instability is induced by the circulatory force represented by the cross-coupled stiffness. The 
nature of this force for this system is expressed by the relationship q = 2 £p<x The value of a 
determines the frequency of the limit cycle in the homogeneous case. Since the fundamental 
driving mechanism for the instability is the same even when excitation is present, it is expected 
that the frequency of a limit cycle under these conditions would remain close to that predicted for 
the homogeneous case. The amplitudes predicted for the homogeneous case would most likely be 
significantly changed by the addition of the unbalance and side-force excitations. However, the 
general trend of increasing amplitude as the global stability limit is approached would be 
expected to hold. 

The unbalance response equilibrium solution demonstrated certain regions in which a 
synchronous sinusoidal response could not exist. In these regions, a combination of synchronous 
response and either limit-cycle instability or subharmonic response might be expected. The 
addition of a side-force excitation would alter the specific range of occurrence and amplitudes of 
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this behavior; however, for some values of side force this behavior might still be expected to 
occur. The stabilizing capacity of the side force has already been discussed. The side force has an 
additional effect of making the system behave like a linear system. This is true for relatively large 
values of side force (compared to deadband and unbalance). 

The subharmonic response analysis already includes the combination of unbalance 
excitation and side-force excitation. However, the form of solution assumed in the procedure 
does not allow for the occurrence of limit-cycle instability. The ranges of possible subharmonic 
response should be valid; however, the existence of this form of solution is not guaranteed. 
Limit-cycle instability may also occur in regions where subharmonic response is possible. The 
results still show the effects that various model parameters will have on the subharmonic if it 
occurs. This is typical of certain types of nonlinear systems where multiple solutions are 
possible. In fact, for a given set of parameters, the initial conditions may determine which 
solution is obtained. Other perturbations, such as the random noise, may also play a strong part 
in the determination. 

The results presented in this section provide general characteristics of the responses of 
the nonlinear rotor system. They provide insights into the effects that the various parameters will 
have and give direction for the simulation studies to be presented in the next section. The various 
behaviors postulated in this section for the general system will be investigated using 
simulations. 


V. SIMULATION RESULTS FOR SIMPLIFIED MODEL 


Numerical results were presented in the previous section for the analytical expressions 
developed in section III. These results consisted of equilibrium response amplitudes, stability 
conditions, and subharmonic response amplitudes. The analyses were developed for specific 
excitation cases (homogeneous, unbalance, and side force) and a specific assumed form of 
solution (subharmonic response). None of the analyses are fully applicable for a system under 
genera] excitation and one whose solution form is not known a priori to be a superposition of 
subharmonics. However, some general conclusions were drawn for the system in the previous 
section for both the restricted cases for which the analyses apply and for the general case. In this 
section, simulation will be used to demonstrate the validity of the results presented and the 
conclusions drawn in the previous section. In addition, results which can only be determined 
through simulation will be presented. These results were obtained using the general turbopump 
model developed in appendix A. The numerical integration method used for the simulation 
solution is discussed in this appendix. 


A. Demonstrations of Equilibria and Stability of Restricted Cases 

The first result developed was the stability condition for the zero deadband case. For 
nominal parameters, this condition was shown to be p < 2.08. This is illustrated (fig. 18) by 
slowly ramping the simulation through this value and observing the divergent growth of the 
response beyond this value of p. A very low level of random noise was used to initiate the 
instability. Fast Fourier transform (FFT) analysis of the response (fig. 19) shows that the 
frequency of the instability is equal to the normalized natural frequency of the system 
(1 radian/second). 
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Figure 18. Demonstration of stability threshold for zero deadband linear model. 
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Figure 19. Magnitude spectrum of y response of figure 18. Spectrum taken at p « 2.1 






The limit-cycle equilibrium radius given by equation (128) is shown in figure 20 for 
nominal values of all parameters. Simulation results are shown on this plot for certain values of 
p. FFT analyses for the various values of p show that the frequency of the limit cycle is 0.48p as 
predicted by equation (121). These spectra are shown in figure 21. 



Figure 20. Homogeneous equilibrium solution amplitude as p is varied. All parameters are 
nominal. Solid line represents analytical solution, circles represent simulation results. 



Figure 21. Cascade spectral plot of simulation of homogeneous equilibrium. 
Spectra taken at dwells in p profile. 
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The mass unbalance equilibrium radius given by equations (78) and (79) is shown in 
figure 22 for nominal values of all parameters. Simulation results are shown for a ramp up to 
p = p/ c . Above this value, the mass unbalance equilibrium results are not valid since ro < 1 and 
the solution is unstable. Simulation results are also presented for the ramp down from this speed, 
illustrating the multivalued solution which exists in a small range near the resonance. In order to 
explore the stability threshold of the unbalance equilibrium (p ue ), a larger unbalance case was 
simulated. This was necessary to generate an equilibrium solution whose value exceeds 1.0. 
Simulation results are presented in figure 23 for the model operating at a steady speed (p = 1.7) 
just below the stability threshold ( p ue = 1.72) and then ramping to and holding at a speed 
(p = 1.73) just above the threshold. A very low level of noise excitation (rj = 0.0001) was used 
to perturb the equilibrium. At the lower speed, the equilibrium is maintained. At the higher speed, 
the amplitude appears to diverge and then limits at a higher level than the equilibrium. FFT 
analysis of the response (fig. 24) shows that the limit cycle instability has emerged along with 
the mass unbalance response. This result shows that at speeds above this stability threshold 
the unbalance response equilibrium cannot be maintained, and a combination of unbalance 
response and limit-cycle instability results; however, it does not show the converse, i.e., it does 
not show that below this threshold the combination response cannot be maintained. It only 
shows that below the threshold the unbalance equilibrium without the limit-cycle instability is 
possible. 



Figure 22. Mass unbalance equilibrium solution amplitude as p is varied. All parameters 
nominal. Solid line represents simulation solution, X symbols represent analytical results. 

The stability analysis of the side-force equilibrium showed that this equilibrium was 
stable for speeds beyond the global stability threshold. This is demonstrated for a nominal case 
(with no unbalance) in figure 25. This case has a low level of noise excitation (r] = 0.001) to 
perturb the system. The predicted stability threshold for this case (equation (116)) is p = 3.0. 
The system clearly remains stable for speeds below this threshold and diverges beyond it. Since 
the analysis only addressed stability in the small, it is of interest to determine how sensitive the 
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Figure 25. Demonstration of side-force equilibrium stability threshold. n = 0 . 001 , a = 0.0, 

gy = 1 .0. Predicted stability threshold is p = 3. 

stability is to perturbation. This was examined by simulating the system at a constant speed of 
p = 2.7 and increasing the amplitude of the noise perturbation (fj ) . Results from this simulation 
are shown in figure 26 and show that a value of rj ~ 0.47 was required to perturb the system 
beyond the range of stability. This is not an exact value due to the random nature of the noise 
excitation and the fact that the value was steadily increasing. An even more interesting result is 
shown in figure 27. This figure presents two cases, one in which the noise is only applied along 
the y axis and the other only along the z axis. The amplitude is increased as before. For the y- 
axis perturbation, the system remains stable for all values of rj y up to 1.0. However, the z-axis 
perturbation behaves just like the dual-axis perturbation. The system becomes unstable above 
ri z ~ 0.4. In hindsight, this should not be a surprising result. The fundamental driving force for the 
instability is a tangential force. The side-force equilibrium point is primarily along the y axis and, 
hence, a z-axis perturbation would impart a tangential velocity to the rotor. One might 
hypothesize then that an unbalance excitation superimposed on the side-force excitation might 
rather easily perturb the system beyond its range of stability. This might be expected since the 
unbalance provides a large, regularly occurring tangential perturbation to the rotor. This will be 
explored in a later section. 

The analytical results presented in section IV showed that the stability threshold for the 
side-force equilibrium increases as the magnitude of the side force decreases. This is true as 
long as the side force is sufficient to cause the displacement to exceed the deadband (ro > 1.0). 
This can be understood by realizing that the greatest asymmetry occurs in the linearized 
stiffness coefficients when the deflection is the smallest. Gunter, 12 among others, has shown the 
stabilizing capacity of asymmetry. However, due to the smaller magnitude of the side force and 


41 







the deflection for this case, this highly stable (in the small) equilibrium might be expected to be 
more sensitive to perturbation than a more highly loaded case. This is shown to be the case in 
figure 28. The result is for a case of side force g y = 0.5 instead of 1.0 as was the case in figure 26. 
The side-force equilibrium stability threshold for this case is p~ 4.25. All other conditions are 
the same and the system becomes unstable when t] = 0.25. For increasing side force, although 
the sensitivity of the stability to disturbance will decrease somewhat, the stability threshold will 
also decrease and, in the limit, will approach the global stability threshold (p g ,). 



ETA BAR 

Figure 28. Sensitivity of side-force equilibrium stability to perturbation, 
p = 2.7, g y = 0.5, r\ y — T]z~ 

Subharmonic response harmonic balance results were presented in section IV for nominal 
values of the model parameters and several variations of each parameter. These results give the 
amplitudes of a subharmonic if it occurs, but they do not address any perturbations or initial 
conditions which may be required to initiate the response. In order to more easily obtain the 
subharmonic response in a simulation, a case with a damping of £= 0.01 will be demonstrated. In 
addition, the destabilizing force parameter cr will be set to zero to avoid potential interaction 
between the subharmonic response and limit cycle instability. The harmonic balance results are 
shown for this case in figure 29. The z-axis component has a much larger amplitude than the y- 
axis component for this case. The maximum value of p for which convergence was attained is 
1.775. The simulation results are shown in figure 30. The transient data in figure 30b has been 
bandpass filtered between the normalized frequencies 0.3 and 1.26 in order to illustrate the 
amplitude of the subharmonic component. The response matches the prediction quite well up to 
p = 1.775. The frequency is exactly one half the excitation frequency, and the amplitudes match 
for both axes. As the speed continues to increase beyond this point, the response begins a 
transition phase where it appears to be seeking a new equilibrium. A new equilibrium is then 
achieved, one for which the harmonic balance procedure failed to converge. Simulation of the 
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same case with noise perturbation added (tj = 0.05) yields a very interesting result (fig. 31). The 
response below p = 1.775 is similar to the previous case. Above this value, however, when the 
response enters the transition phase it does not arrive at the same equilibrium that it did 
previously. Instead, a harmonic equilibrium (i.e., no subharmonic response) is obtained. This 
result suggests that the subharmonic solution above p = 1.775 is not very stable. This might 
explain why the harmonic balance procedure fails to converge to a solution in this speed range. 



Figure 29. Subharmonic response solution versus p. <r = 0.0, £= 0.01. 

Magnitudes of y- and z-axis components. 

B. Interaction Between Limit Cycle and Mass Unbalance Response 

The previous section dealt with the restricted cases for which equilibria and stability have 
been determined analytically. These cases assume a specific excitation form and/or a specific 

form for the solution. This section and those that follow will deal with cases for which the 
assumptions do not apply. 

The mass unbalance equilibrium simulation was presented in figure 22. This case was 
restricted to speeds below p, c since this equilibrium (r 0 > 1.0) would be unstable beyond that 
point. Proceeding beyond this point results in a combination of unbalance response and limit cycle 
instability. This is shown for the same model parameters in figure 32. The initiation of the limit 
cycle is most clearly evident from the cascade spectral plot. The frequency of the limit cycle is 
^ approximately equal to <xp; however, the unbalance excitation does alter it somewhat. 
Modulation frequencies of the excitation frequency and the limit cycle frequency can also be 
observed in the spectral data. As predicted by Day, 14 these frequencies can occur at all multiples 
of the difference between the two, plus or minus the limit cycle frequency. The system in this 
case possesses a unique solution. For p < p lc , the only possible solution is the unbalance 
equilibrium. For p > p/ c , the only possible solution is the combination solution. 
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Another mass unbalance equilibrium was presented in figure 23. This solution is for a 
case of increased unbalance (a = 1.0). For this case, the equilibrium radius exceeds the deadband 
(ro > 1.0). This allows for a stable equilibrium beyond the previous limit of p = pi c . The stable 
equilibrium is possible for speeds up to the equilibrium stability limit p ue = 1.72. This was 
demonstrated in figure 23. Although the stable equilibrium is possible up to p = p ue , it is not 
guaranteed to be a unique solution. For values of p greater than p/ c , a combination of unbalance 
response and limit cycle instability is possible even though a stable equilibrium is also possible. 
This is illustrated in figure 33 by beginning where the simulation of figure 23 ended and ramping 
back down into the potentially stable region. The limit cycle is maintained in this region down to 
p ~ 1.45. To verify that the response is not merely the transient decay of an unstable response at 
higher speeds, the ramp down is stopped at p = 1.5 and the limit cycle continues (fig. 34). The 
absence of the limit cycle below p = 1.45 suggests that the unbalance excitation introduces a 
threshold between p/ c and p ue below which the combination response is not possible. The limit 
cycle can also be obtained in this region by perturbation. Figure 35 shows a case where random 
noise excitation is used to perturb the system (rj= 0.1). The limit cycle response is initiated in 
this case at about p = 1.6 which is below p ue . As speed continues to increase to p = 2.0, the limit 
cycle instability increasingly dominates the response. The amplitude grows as predicted by the 
homogeneous analysis and the frequency becomes equal to op. 
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Figure 33. Nonuniqueness of unbalance equilibrium solution. Ramp down from 
p = 1.73. pic = 1.04, p ue = 1.72, n = 0.0001. 
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Figure 34. Repeat of simulation of figure 33 with dwell on ramp down at p = 1.5 

(dimensionless time equal 24,150). 


C. Interaction Between Limit Cycle and Side-Force Response 

*-ie stabilizing capacity of the side force was explored in an earlier section. The effects of 
side force on limit cycle instability will be addressed here. Since the side force can stabilize the 
system at speeds greater than the global stability threshold, one might speculate that it could 
also suppress the limit cycle instability that can occur at speeds below this threshold. Examples 
of this are presented in figures 36 through 40. For each of these figures, the simulation has no 
unbalance excitation. The side force is initially zero and is slowly ramped up to a maximum and 
then back down to zero at the same rate. 

Figures 36 through 38 represent three different values of p: 1.25, 1.75, and 2.0, 
respectively. The maximum side force for each is 3.0. The simulation was initiated by quickly 
ramping to the operating speed and allowing the limit cycle to achieve steady-state conditions 
before the side force was applied. For the first two cases, the response exhibited a hysteretic 
behavior. While increasing the side force, a value was reached which caused the limit cycle to 
cease. While decreasing the side force, the limit cycle remained suppressed until a lower value of 
side force was reached. The more unstable case (p = 1.75) required more side force to suppress 
the limit cycle on the up ramp than did the case with p = 1.25. For the third case (p = 2.0), the 
maximum side force (g y = 3.0) was not sufficient to suppress the limit cycle. This case was 
suppressed for a side force of approximately 6.0. With the larger maximum side force, the system 
exhibited the same hysteretic behavior as in the other cases. It is interesting to note that the 
limit cycle reinitiates at about the same value of side force on the down ramp for all three cases. 
This corresponds to the value for which the side force is insufficient to displace the rotor beyond 
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Figure 36. Limit cycle suppression by side force. Maximum g y = 3.0, 
p= 1.25,a = 0.0, r? =0.0001. 
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Figure 37. Limit cycle suppression by side force. Maximum g y = 3.0, 
p= 1.75,0 = 0.0, n =0.0001. 
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Figure 40. Effects of noise perturbation on limit cycle suppression by side force. 
Maximum g y = 3.0, p = 1.75, a = 0.0, 0 = 0.1, and rj = 0.5. 


the deadband. Although the limit cycle suppression occurs at different values of side force, in 
each case it occurs approximately at the value where the rotor trajectory passes within the 
deadband zone. In other words, if the magnitude of the deflection vector falls below 1.0 during a 
portion of each period of the limit cycle, suppression of the limit cycle is imminent. The sensitivity 
of the suppression and reinitiation thresholds to noise perturbation is also of interest. Figure 40 
shows the same system as figure 37 for two cases of noise excitation, Tj = 0.1 and q = 0.5. The 
first case shows little effect from the noise. The second shows a small reduction in the side force 
required to suppress the limit cycle (from =1.15 to =0.95). It also shows a small increase in the 
reinitiation threshold (from =0.3 to =0.4). One final observation that should be made for these 
cases is that the frequency of the limit cycle is relatively unaffected by the magnitude of the side 
force. This can be seen from the cascade spectral plot of the p = 1.75 case shown in figure 41. 


D. Interaction Between Limit Cycle and Subharmonic Response 

The capacity of mass unbalance and side-force excitation to inhibit limit cycle instability 
has been explored for each excitation individually. The effect of these excitations applied 
simultaneously will now be explored. One effect that is anticipated is the entrainment of the limit 
cycle frequency by the subharmonic response frequency. This is only possible for this system 
when both excitations are present since both are required to produce the subharmonic response 
phenomenon. 

The approach now taken is to repeat the previous numerical experiments (figs. 36 through 
39) with mass unbalance added. This will be done for two values of the unbalance parameter, a = 
0.5 and a = 1.0. The effects of noise will also be explored. Figures 42 through 44 present the 
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Figure 41. Cascade spectral plot of response of figure 37. Spectra taken in dimensionless 

time increments of 400. p = 1.75. 
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Figure 42. Interaction between limit cycle and subharmonic response. 
Maximum g y = 3.0, p = 1.25, a = 0.5, and rj = 0.0001. 
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Figure 44. Interaction between limit cycle and subharmonic response. 
Maximum g y = 10.0, p = 2.0, a = 0.5, and H = 0.0001. 








results for a = 0.5. In figures 42 and 43, p has values of 1.25 and 1.75, respectively, and g y has a 
maximum value of 3.0 as before. In figure 44, p has a value of 2.0 and g y has a maximum value of 
10.0. The results for these cases are similar to the side force only cases. One major difference is 
the entrainment of the limit cycle by the subharmonic response. In figure 42, this occurs almost 
instantly as the side force is applied. Comparing figure 43 with the corresponding case with no 
unbalance (fig. 37), it is clear that the limit cycle transitions to subharmonic at a value of g y which 
is less than the value which suppressed the limit cycle in figure 37 (=0.75 versus =1.25). In 
addition, the subharmonic response is maintained beyond the limit cycle suppression value (up to 
g y ~ 1.75). Upon decreasing the side force, the subharmonic and limit cycle remain suppressed 
until a lower value of side force is reached. Adding noise affects this behavior as shown in figure 
45 (*7 = 0.5). The reinitiation occurs for larger values of side force than without the noise. The 
cascade spectral plot for this case (fig. 46) shows the distinct frequency shift that occurs when 
the limit cycle becomes entrained. The case with p = 2.0 (fig. 44) does not demonstrate the 
entrainment. This case behaves almost identically to the corresponding case without unbalance 
(fig. 39) with a small harmonic component superimposed due to the unbalance excitation. The 
absence of the subharmonic is due to the fact that p = 2.0 is above the range of possible 
existence of the subharmonic, as shown in figure 13. 


The results for the same cases with increased unbalance ( a = 1.0) exhibit somewhat 
different behavior. The low speed case (p = 1 .25) does not exhibit any limit cycle instability. This 
is due to the existence of a stable unbalance equilibrium for this speed. The lack of a one-half 
subharmonic is due to the fact that p = 1.25 is below the range of possible existence. As the side 
force increases, a low amplitude two-thirds subharmonic develops for a small range of side force 
when noise excitation is present (ff = 0.5). The results are shown in the cascade spectral plot in 
figure 47. The intermediate speed case (p = 1.75) exhibits subharmonic response almost 
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Figure 45. Interaction between limit cycle and subharmonic response. 
Maximum g y = 3.0, p = 1.75, a = 0.5, and rj = 0.5. 
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Figure 46. Cascade spectral plot of response of figure 45. Spectra taken in dimensionless 
time increments of 400. p 1.75, a = 0.5, and f? = 0.5. 
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Figure 47. Interaction between limit cycle and subharmonic response. Cascade spectral plot of 
response. Spectra taken in dimensionless time increments of 400. Maximum g y = 3.0, 

p = 1.25, a = 1.0, and rj = 0.5. 






immediately after the initiation of the side-force ramp. Due to unbalance excitation effects, the 
limit cycle frequency ratio is very near 0.5 for this case. The subharmonic is suppressed when the 
side force exceeds the value for which existence is possible. The subharmonic reinitiates at a 
lower value when the side force is decreased. These results are shown in figure 48. The high 
speed case (p = 2.0) again exhibits only limit cycle instability since it falls beyond the range of 
possible subharmonic. The results are very similar to the previous cases ( a = 0.0 and a = 0.5) 
and are shown in figure 49. 



Figure 48. Interaction between limit cycle and subharmonic response. Cascade spectral plot of 
response. Spectra taken in dimensionless time increments of 400. Maximum g y = 3.0, 

p = 1.75, a = 1.0, and V = 0.5. 
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Figure 49. Interaction between limit cycle and subharmonic response. Cascade spectral plot of 
response. Spectra taken in dimensionless time increments of 400. Maximum g y = 10.0, 

p = 2.00, a = 1.0, and rj = 0.5. 
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E. Effects of Operating Profile on Interacting Responses 

The effects of the unbalance, side force, and combined unbalance and side force on the 
limit cycle instability have been demonstrated. The results have shown a strong dependence on 
initial conditions and perturbations. This is indicated by the hysteretic behavior during the side- 
force ramp up and back down, for example, and by the sensitivity to noise perturbation. These 
results were obtained using certain numerical experiments where parameters were intentionally 
varied in order to gain insight into the response of the system. These variations were not 
necessarily representative of actual operating profiles of a turbomachine. Due to the sensitivity of 
the system to initial conditions, it is important to consider realistic operating profiles to assess 
the potential for the occurrence of limit cycle behavior or subharmonic response. This will be 
accomplished using a series of simulations where speed is ramped up and back down while the 
other system parameters are held fixed. The simulations are initialized by rapidly ramping speed 
from zero to a value of p — 1.0. Speed is then held constant for a period of time. The side force 
begins at zero initially and ramps to its steady value somewhat slower than speed. It reaches its 
maximum at the end of the steady hold time for speed. Speed is then ramped up to p = 2.0 and 
back to p = 1.0 while the side force remains constant. These profiles are shown in figure 50. This 
approach allows the system to achieve a steady-state initial condition at a speed which is 
outside the range of limit cycle or subharmonic response. The maximum value of the side force, 
the unbalance, and the random noise parameter are the model parameters that will be varied in 
this series of experiments. 

The first three cases examined all have unbalance magnitudes of a = 0.5. The first has a 
side force magnitude of g y = 1.5. For the second g y - 1.0 and for the third g y = 0.5. Results for the 
first case are shown in figure 51. This case also has random noise perturbation (rj = 0.5). No 
limit cycle or subharmonic occurs for this case. A low-level response of the system resonance to 
the noise excitation is visible in the plot. Without the noise excitation this case exhibited only 
synchronous harmonic response to the unbalance excitation. This case represents a case similar 
to that of figures 43 and 45. In the earlier simulations, speed was held fixed and the side force 
was varied while the converse is true in the current simulation. Examining the response in figure 
45 at a value of g y = 1.5 shows that a subharmonic existed on the up ramp but only the 
synchronous harmonic existed at this value on the down ramp. These results indicate that the 
initial conditions and/or the perturbation used in the current simulation were not suitable for 
initiating the subharmonic response. The second case (g y = 1.0) also failed to exhibit any 
nonsynchronous response in the absence of noise. With noise excitation (f) = 0.5) however, both 
subharmonic response and limit cycle instability are observed. These results are shown in figures 
52 and 53. The subharmonic initiates at p = 1.6. The response transitions to limit cycle instability 
above p= 1.8. This is approximately the upper limit of possible subharmonic response for this 
case. For the case of low damping and no cross-coupled stiffness examined previously (figs. 29 
and 30), the response in this region appeared to be seeking a subharmonic equilibrium that was 
not predicted by the harmonic balance procedure. When noise excitation was added (fig. 31) this 
did not occur. The current results suggest that for a system which has cross-coupled stiffness of 
the type addressed here (cr near 0.5), the response will transition to the limit cycle instability in 
this region and not arrive at the weakly stable subharmonic solution. This would be expected 
since the system is moving very close to the global onset speed of instability p g i. The third case 
(£ = 0.5) behaves very similarly to the second. The results are shown in figures 54 and 55. One 
significant difference is that noise excitation was not necessary to initiate the subharmonic 
response. The lower value of side force causes the displacement to fall within the deadband 
(r < 1.0) which perturbs the system to such an extent that the subharmonic is initiated. Once the 
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subharmonic is initiated, it transitions to the limit cycle when p exceeds the range of possible 
subharmonic. The only other difference noted is in the specific values of p at which transition 
occurred. 

The three cases just examined were repeated with an unbalance value of a = 1.0 instead 
of a = 0.5. The results are shown in figures 56 through 61. For each of these cases, the 
subharmonic was initiated and transitioned to limit cycle instability without the noise 
perturbation. When the side force is increased to = 2.0, noise is necessary to cause the 
subharmonic response to initiate (figs. 62 and 63). 
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Figure 56. Effects of operating profile on interacting responses. g y = 1.5, a = 1.0, and H = 0.0001. 



0-00 0-25 0-S0 0.75 1-00 1.25 I SO 1.75 2.00 2-25 2-50 


NORMALIZED FREOUENCr 

Figure 57. Cascade spectral plot of figure 56. Spectra taken in dimensionless time increments of 

400. g y = 1.5, a = 1.0, and t? = 0.0001. 
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F. Effects of Other System Parameters on Interacting Responses 

The effects of speed, side force, mass unbalance, and random noise on the interacting limit 
cycle and subharmonic response have been explored. The effects of the cross-coupled stiffness 
parameter cr, the nonlinear stiffness parameter y, and the damping parameter £ are also of 
interest. The previous case with g y = 1.0 and a = 1.0 (figs. 58 and 59) will be used as the basis 
for making variations in these other system parameters. 

The effects of a were examined by changing its value from the nominal of 0.48 to 0.40 and 
0.52. Results of the first case are shown in figures 64 and 65. The maximum speed was increased 
to p = 2.5 for this case since that is the value of the global instability threshold. Random noise 
was used to increase the probability of initiation of potential subharmonic response or limit cycle. 
A subharmonic response occurs for this case in the speed range where it can exist. When this 
range is exceeded, the subharmonic disappears and no limit cycle appears. This is in contrast to 
the previous cases where the subharmonic transitioned to a limit cycle. For this value of cr 
however, the destabilizing force is insufficient to overcome the stabilizing capacity of the side 
force. As speed is increased and the system becomes less stable, the limit cycle initiates at a 
frequency of ~o p. This occurs at a higher speed than in the cases with the nominal value of <7. 
There are, therefore, two distinct regions of speed where subharmonic response and limit cycle 
instability individually occur when o= 0.40. The case with <r = 0.52 is shown in figures 66 and 67. 
The behavior here is similar to the baseline case except for the obvious difference in the limit 
cycle frequency (~0.52p). One other ditference is that the global instability limit (Pgi — 1*92) is 
just slightly above the upper limit of possible subharmonic response. This causes the transition 
from subharmonic to limit cycle to occur just before the onset of divergence. This simulation was 
executed with a ramp rate half that of the previous cases in order to identify the transition. 
Another difference is that the transition back to subharmonic response on the ramp down exhibits 
considerably more hysteresis than the previous cases. In other words, the range of speed for 
which either response can occur is greater. 

The effects of /were examined by changing its value from the nominal of 0.75 to 0.50 and 
0.25. These values represent progressively more linear cases. The results for y= 0.50 are shown 
in figures 68 and 69. The general characteristics of the response are the same as the nominal 
case. As predicted by the subharmonic response analysis, the subharmonic remains for higher 
speeds and has a lower amplitude. As predicted by the homogeneous limit cycle analysis, the 
limit cycle has a lower amplitude also. The results for y= 0.25 are shown in figures 70 and 71. 
This case did not develop either subharmonic response or limit cycle. This is due to the very 
narrow range of speed and low amplitude for potential subharmonic and the low amplitude for 
potential limit cycle. Since this system is nearly linear, this might also have been intuitively 
expected. 

The effects of £ were examined by changing its value from the nominal of 0.10 to 0.05 and 
0.20. It should be noted that when £ is changed the cross-coupled stiffness is also changed 
proportionately according to equation (125). The results for £ = 0.05 are shown in figures 72 and 
73. The only difference noted is that the transition from subharmonic to limit cycle occurs at a 
higher speed than in the nominal case. This is due to the fact that the upper limit of potential 
subharmonic is higher for this case (fig. 16). The results for £ = 0.20 are shown in figures 74 and 
75. As might be expected, the transition occurs at a lower value of speed for this case. 
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Figure 66. Effects of other system parameters on interacting responses. 
c — 0.52, g y = 1.0, a = 1.0, and Tj = 0.5. 
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Figure 67. Cascade spectral plot of figure 66. Spectra taken in dimensionless time increments 

of 400. <r= 0.52, g y = 1.0, a = 1.0, and rj =0.5. 
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Figure 68. Effects of other system parameters on interacting responses. 
Y= 0.5, g y = 1.0, a - 1.0, and rj = 0.5. 
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Figure 69. Cascade spectral plot of figure 68. Spectra taken in dimensionless time increments 

of 400. y= 0.5, g y = 1.0, a = 1.0, and fj = 0.5. 


1 








Figure 70. Effects of other system parameters on interacting responses. 
y= 0.25, g y = 1.0, a = 1.0, and rj = 0.5. 
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Figure 71. Cascade spectral plot of figure 70. Spectra taken in dimensionless time increments 

of 400. 7 = 0.25, g y = 1.0, a = 1.0, and r\ = 0.5. 
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Figure 74. Effects of other system parameters on interacting responses. 
£ = 0.20, g y = 1.0, a = 1.0, and tj =0.5. 





Figure 75. Cascade spectral plot of figure 74. Spectra taken in dimensionless time increments 

of 400. £=0.20, g y = 1.0, a = 1.0, and H =0.5. 
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G. Nonuniqueness of Solutions 


The nonuniqueness of either the subharmonic response solution or the limit cycle solution 
when mass unbalance and side force excitations are applied simultaneously has already been 
indicated by the hysteretic behavior of some of the simulations. This occurred for both side-force 
ramps and speed ramps. Results for another simulation which clearly illustrates this are shown 
in figures 76 and 77. This case is the same as that shown in figures 58 and 59 (a = 1.0, g y = 1.0, 
£, a, and /nominal) with two exceptions. First, random noise (rj = 0.5) has been added. 
Second, the speed profile has been altered to dwell on the up ramp at p = 1.845 and remain at 




Figure 76. Nonuniqueness of interacting response solution. Dwell in operating profile at 
p = 1.845 (dimensionless time equal 12,000). g y = 1.0, a = 1.0, and rj = 0.5. 
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Figure 77. Cascade spectral plot of figure 76. Spectra taken in dimensionless time increments 

of 400. g y = 1 .0, a = 1 .0, and 7? = 0.5. 
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that value for the remainder of the simulation. From the Figures it can be seen that the response 
is initially a subharmonic. After some time it changes to limit cycle. Some time after this it 
changes back to a subharmonic. All model parameters are fixed during this time penod. The 
random noise has the effect of superimposing a variation on the side force. This behavior o 
changing back and forth between solutions when no apparent change has occurred in the system 
is very significant when reviewing test data. This will become more apparent when test data 
from the SSME HPFTP is presented in section VI. 


H. Effects of Mass Unbalance on the Stability of Side-Force Equilibria 

The stability of the side-force equilibrium was analyzed in section III and demonstrated 
with simulations earlier in this section. The possible effects of mass unbalance perturbing t e 
stability of this equilibrium were discussed. The effects have been explored using simulations of 
the nominal model with g y = 1.5 and two values of unbalance ( a - 0.5 and a - 1.0). The stabi y 
threshold for the side-force equilibrium in this case is p - 2.7. Speed was ramped from p 1.0 
p = 3.0 after initiating the simulations as discussed previously. Results for both cases are 
presented in figure 78. Results for a = 0.5 show that this value of unbalance does not perturb the 
equilibrium beyond its range of stability. It remains stable for values of p up to 2.7. For a - I.O, 
the unbalance is sufficient to initiate a subharmonic response at p- 1.6. As speed increases, this 
transitions to limit cycle instability at p = 1.9. Now the system is in a limit cycle instability when 
the speed increases beyond the global instability threshold (p gi - 2.08) and the response 
diverges. The effect of the unbalance is to alter the state of the system when the global 
instability threshold is crossed. In another sense, the subharmonic response causes the initiation 
of the limit cycle which diverges when the threshold is crossed. If the unbalance and side-force 
parameters are such that the subharmonic does not occur, the side force can stabilize the system 
beyond the global threshold even in the presence of unbalance and noise perturbation. 
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Figure 78. Effects of mass unbalance on the stability of side-force equilibria. Stability threshold 

is p = 2.7. g y = 1.5, rj = 0.0001, a = 0.5, and a = 1.0. 
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The behavior of the model under general loading conditions has been characterized in 
terms of its behavior under certain restricted conditions. These restricted conditions are the 
homogeneous response (which describes the limit cycle behavior), the mass unbalance 
equilibrium, the side-force equilibrium, and the general loading case with the assumption of 
harmonic and subharmonic response. These results will be extended to the SSME HPFTP in 
section VI. 


VI. EXTENSION OF RESULTS TO COMPLEX MODEL 


The previous sections have dealt with a simplified model of a turbopump in order to 
enhance the understanding of the basic phenomena being studied. The objective of this section is 
to demonstrate that these results extend to a more complex, realistic model of an actual 
turbopump. This will be achieved by examining the HPFTP of the SSME. Test data will be 
presented for cases where the phenomena appear to have occurred. Analyses and simulations of 
a model of the HPFTP will also be discussed. 


A. Description of SSME HPFTP 

The SSME is manufactured by Rockwell International, Rocketdyne Division for NASA. It 
is a liquid hydrogen/liquid oxygen staged combustion rocket engine. The primary components of 
such an engine are the main combustion chamber and nozzle, the high pressure turbopumps 
which feed the fuel and oxidizer to the main combustion chamber, and the combustion system 
which drives the turbopumps. In the case of the SSME, each pump has its own combustion 
device known as a prebumer. The fuel is partially burned in these prebumers (oxygen/hydrogen 
mixture ratio approximately one to one) and the resulting combustion gases drive the turbines. 
These gases then proceed to the main chamber where they are completely burned according to 
stoichiometric balance (mixture ratio approximately six to one). A schematic diagram of the 
propellant flow is shown in figure 79a. The HPFTP is shown in figure 79b. For engine operation 
at 109 percent of rated power level, the HPFTP runs at approximately 36,600 r/min, depending on 
actual engine performance. It consists of a three-stage centrifugal pump section which is driven 
by a two-stage turbine. The rotor is supported primarily by two pairs of angular contact ball 
bearings, one pair on each end. The two-pump interstage seals also provide significant restoring 
forces for relative lateral rotor motion. These seals provide the majority of the damping and 
cross-coupled stiffness forces. The turbine section provides significant additional cross-coupled 
stiffness due primarily to the Alford effect. The Alford effect is a variation in the aerodynamic 
efficiency of individual turbine blades as the turbine disk moves eccentric with respect to the 
stator. The blades on one side will be more efficient than those on the other due to their smaller 
tip clearance. The net effect is a tangential force proportional to the radial displacement which is 
modeled as a cross-coupled stiffness. Additional stiffness, damping, and cross-coupled stiffness 
forces arise at the turbine interstage seal. The impellers and their associated seals produce 
forces as well; however, their magnitudes are much smaller than the pump interstage seal forces 
These forces will be neglected here in order to simplify the model. The error introduced by doing 
so is no greater than the error due to the uncertainty in the dominant pump interstage seal and 
ball bearing forces. 

The nominal data used to define the HPFTP model used in this study are contained in 
appendix B. These data were adjusted within their uncertainties in order to achieve the desired 
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behaviors as observed in tests. The adjustments are discussed in the section dealing with model 
results. These data were provided by the turbopump manufacturer (Rocketdyne) with the 
exception of the pump interstage seal cross-coupled stiffness. This data has been replaced by a 
function of the form of equation (125). The damping value supplied by the manufacturer is 
multiplied by the shaft angular velocity and then by the destabilizing force parameter cr. This 
simplifies the process of relating the HPFTP model results to the simplified model results. 

The model for the HPFTP is more complex than the single mass model in several 
respects. One of the most obvious differences is the addition of the housing dynamics to the 
model. Another distinction is that all parameters are distributed along the axis of the rotor. The 
support characteristics (both linear and nonlinear), rotor and housing mass and stiffness 
properties, and the excitation forces all have independent (discrete) distributions along the rotor, 
in general. These differences are complicating enough, however, the most significant complexity 
is not as obvious. Virtually every parameter in the model varies with turbopump operating speed 
either directly or indirectly. The rotor and housing free-free dynamic characteristics and the mass 
unbalance distribution are exceptions. The stiffness and damping coefficients and the side-force 
excitation vary either with speed directly or with engine power level (which can be related to 
speed). The clearances, geometric eccentricities from centerline, and random noise excitations 
would be expected to vary in the actual machine. These variations are unknown, however, and 
are not prescribed in the model. The fact that they probably occur must be recalled when 
interpreting test data and simulation results. Finally, simply the number of parameters in the 
model make interpretation of results difficult. The same result or trend could probably be obtained 
with more than one set of parameters (or "recipe"). 


B. HPFTP Test Data 

Development of the SSME began in 1971. High power level testing began in 1978. Since 
that time hundreds of engine tests have been performed yielding massive amounts of data. 
During the course of the engine’s development, the design of the HPFTP has evolved and there 
are many different configurations that have been tested. The subsynchronous vibration problems 
that motivated this work have changed with the design changes. An in-depth review of the 
history of the problems and the design changes would take volumes, however, Hawkins 17 
provides a good summary. There are three primary conclusions with relevance to this work which 
can be drawn from the historical data. First, the occurrence of the subsynchronous vibration, its 
amplitude, and its frequency are erratic. Only a certain percentage of all tests exhibit the 
phenomena and the frequencies and amplitudes vary. The frequencies fall within the range of 47 
to 56 percent of rotor speed for all configurations and 47 to 52 percent for the current 
configuration. Second, configurations with higher side forces do not exhibit the phenomena as 
often and the amplitudes tend to be lower. The side-force differences are due to changes in the 
turbine discharge section which altered the pressure distribution in the turbine section. Third, the 
increased stiffness of the currently used pump interstage seals reduces the amplitudes of the 
phenomena. The initial seals used in high power level testing were three step, smooth seals and 
the current seals are straight, smooth seals. One might expect that these configurations would 
also have a lower frequency of occurrence; however, they were generally tested to higher power 
levels at which they are less stable. 

In most tests, the rotordynamic instrumentation on the HPFTP has been limited to 
accelerometers mounted externally on the housing. This makes it somewhat difficult to determine 
what dynamic behavior of the rotor might be creating a particular vibration response on the 
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housing. Specially instrumented units have been built on occasion which have internal rotor 
displacement measurements. One test of one of these units has been selected for discussion 
here. The HPFTP in this test displayed significant levels of subsynchronous vibration with very 
interesting characteristics. The test was designated as 750-270 and the HPFTP was unit 
number 2708R1. The turbopump was instrumented with radial displacement measurements in 
two axes at the seal between the first and second pump stages. It also contained the normal set 
of external accelerometers. 

The power level profile for test 750-270 is shown in figure 80. After throttling down to 80- 
percent power level, the level was slowly increased to 109-percent power level. The actual 
profile was a series of steps of one-half percent power level with 3-s dwells at each step. Upon 
reaching 109-percent power level at 238 s, the level was held constant until 271 s. During this 
period, the liquid oxygen tank pressure was reduced (vented) to simulate flight conditions. This 
lowers the liquid oxygen inlet pressure to the low pressure oxidizer turbopump and, hence, 
lowers the liquid oxygen outlet pressures throughout the system. The engine controller 
compensates for this perturbation by adjusting valves in the system and maintains system thrust 
(power level) and mixture ratio. The liquid oxygen vent schedule is superimposed on the power 
profile in figure 80. 

The response of the HPFTP to the power profile described above is shown in figures 81a 
through 81g. The figures contain frequency spectra of one of the radial displacement 
measurements. These data were provided by the engine manufacturer (Rocketdyne). Many 
measurements and many methods of processing the data are available. This measurement, is 
representative of the response of the HPFTP, and this method of processing yields a concise 
representation of the important features of the response. This series of plots presents a 
sequence of frequency spectra beginning at 170.8 s and continuing until 272.7 s. Each spectrum 
represents nine spectral averages (to reduce noise). The frequency resolution of the spectra is 
5.0 Hz, therefore, each individual spectrum requires 0.20 s of data for analysis. Each plot of the 
average of nine spectra, then, represents 1.80 s of data. The samples are taken at 3-s intervals 
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Figure 81a. Test 750-270 displacement measurement spectral data. Beginning 170.8 s after engine start. 
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Figure 81b. Test 750-270 displacement measurement spectral data. Beginning 185.8 s after engine 










Figure 81c. Test 750-270 displacement measurement spectral data. Beginning 200.7 s after engine start. 
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Figure 8 Id. Test 750-270 displacement measurement spectral data. Beginning 215.8 s after engine 









Figure 81e. Test 750-270 displacement measurement spectral data. Beginning 230.7 s after engine start. 
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Figure 8 If. Test 750-270 displacement measurement spectral data. Beginning 245.7 s after engine start. 
















in order to correspond to the dwells in the profile. The spectra are arranged in the figures from 
bottom to top for increasing time. The distinct spike near 600 Hz is the synchronous component 
due to mass unbalance excitation. A line has been superimposed on the plots at the frequency 
corresponding to one-half this synchronous frequency. This is to highlight any occurrence of 
subharmonic response. 

The spectra in figure 81a are typical of responses earlier in the profile. The noise in the 
system is exciting the first resonance; however, neither subharmonic response nor limit-cycle 
instability appear to be present. It should be noted that the vertical axis scale is logarithmic for 
these plots. As the power level increases, figures 81b through 81e show the inception of one-half 
order subharmonic response. The subharmonic increases in amplitude throughout this interval. 
The last two plots in figure 81e occur after the power level is held constant at 109 percent. All 
plots in figure 8 If are also at this constant power level. The first plot in this figure shows that the 
system is still exhibiting subharmonic response. The remaining plots indicate a transition 
between subharmonic response and limit-cycle instability. This is very similar to the behavior 
described in section V (figs. 76 and 77) for the simplified model. The amplitude increases 
significantly as the transition occurs. The modulation frequency of the excitation frequency and 
the limit cycle frequency which was discussed in section V (fig. 32) is evident during the 
transition. Figure 81g shows the fully developed limit-cycle behavior. The frequency of the limit 
cycle is approximately 47 percent of the rotation frequency. The last two plots in this figure occur 
after power is reduced from 109 percent and are not of interest here. It should be noted that the 
liquid oxygen venting was initiated during the period of fixed power level when the transition from 
subharmonic response to limit cycle occurred. Although an exact mechanism is not known, the 
venting may have created a small perturbation to the HPFTP which initiated the transition. Also, 
as demonstrated in section V, the system noise itself may have been sufficient to initiate the 
transition and the timing with the vent profile may have been just coincidence. 

The test data described above can be interpreted in the context of the analytical and 
numerical results of the previous sections. The machine appears to be operating at the upper limit 
of the range of possible subharmonic response. For the simplified model, this upper limit ranged 
from about 1.7 to about 1.9 times the zero deadband linear system resonance, depending on the 
system parameters. In addition, the transition to limit cycle instability with higher vibration 
amplitudes indicates that the machine is approaching the global onset speed of instability as 
defined in section III. For the simplified model, the 47-percent limit cycle frequency ratio would 
imply that the global onset speed would be 2.13 times the linear system resonance. For the 109- 
percent power level speed of approximately 36,300 r/min, using the upper limit of possible 
subharmonic given above yields a linear system resonance between 19,105 r/min (318 Hz) and 
21,353 r/min (356 Hz). The corresponding range of onset speeds of instability is 40,694 r/min to 
45,482 r/min. These values will be used as a guide in assessing the results from the 
mathematical model. 


C. Linear Analysis Results 

The mathematical model results can be separated into linear analysis results and 
nonlinear simulation results. The primary linear analysis results of interest here are the 
eigenvalues and eigenvectors of the linear system obtained from the zero deadband assumption. 
Since the parameters of the system model are functions of operating speed, the eigenvalue 
problem must be solved for many speed values in the range of interest. For each different value of 
speed, the entire set of eigenvalues and eigenvectors (in this case, 39 complex conjugate pairs) 


87 



are recalculated. In general, since the system matrix is different at each speed, the eigenvalues 
and vectors for one speed may not be related to those at another. However, there is a known 
relationship between the elements of the system matrix at different speeds. Because of this, one 
would expect to be able to relate the system characteristics for one speed to those at another 
speed which is relatively close to the first. This is in fact the case, however, it is not entirely 
straightforward. As speed is varied, the eigenvalues tend to move along loci whose patterns 
become apparent by visual observation. The associated eigenvector of an eigenvalue at one 
speed for a particular locus is usually very similar to that at another speed on the same locus. 
This is not always the case, however. For large speed changes, the dynamics of the system can 
change dramatically, and there may be no recognizable relationship between eigenvectors at one 
speed and those at another. This may be true even though the eigenvalues traced out clearly 
identifiable loci when migrating as speed varied. Another difficulty is when two loci intersect or 
nearly intersect. In these instances, the characteristics of the eigenvectors associated with the 
loci may switch. This type of behavior is observed in the results from the HPFTP model. When 
this occurs, the only meaningful association between eigenvalues at one speed and those at 
another must be based on similarity of their associated eigenvectors, not on the patterns of their 
loci. 


The eigenanalysis results can be presented in several ways. One way is to plot the 
eigenvalues in the complex plane where the locus patterns are evident (root locus). This 
presentation has the disadvantage that the corresponding speeds are not readily visible. Another 
way is to plot the real and imaginary parts of the eigenvalues versus speed. This presentation 
does not provide as much visibility to the loci patterns but clearly indicates the speed 
correspondence. Yet another method is to plot the critical damping ratio associated with each 
pair of complex roots versus speed. This provides some "calibration" for the real part of the 
eigenvalue indicating its relative stability. The complex eigenvectors are presented in the manner 
described in appendix A. Due to the large number of eigenvalues for this system, eigenvalue and 
eigenvector information will only be provided for those eigenvalues related to the limit cycle 
instability and subharmonic response. 

The linear analysis (zero deadband assumption) results for the HPFTP using the nominal 
data discussed in appendix B are presented in figures 82 through 84. These figures present the 
root loci, imaginary components, and critical damping ratio, respectively, for two eigenvalues 
designated No. 3 and No. 4. On the root loci (fig. 82), the "X" symbol indicates the starting speed 
(10,000 r/min) and each circle corresponds to a 1,000-r/min increment. These eigenvalue loci 
exhibit the switching behavior discussed previously. At the higher speeds, the eigenvector 
associated with eigenvalue No. 3 possesses the characteristics of the unstable behavior. This 
eigenvector is shown in figure 85. This figure displays the relative component (rotor minus 
housing) of the eigenvector. The rotor precesses in a forward direction (except at the pump end 
bearings) with only a small amount of housing motion. The motion can be described as a rigid 
body translation with superimposed rotor flexing. From figures 82 and 84, it can be seen that the 
locus for this eigenvector is moving toward instability. From figure 83, the frequency is observed 
to vary between 250 and 360 Hz with a value of about 335 Hz at the 109-percent power level 
speed (36,300 r/min). At the lower speeds, the eigenvector associated with eigenvalue No. 4 
possesses the characteristics of the unstable behavior. The locus for this eigenvalue heads 
toward the right half plane initially. The critical damping ratio reaches a minimum at 25,000 r/min 
and the locus then reverses direction. This is due to the switching which takes place between 
this locus and the one for eigenvalue No. 3. 


88 



M 

R 

G 


1 

N 

R 

R 

Y 

P 

n 

R 

T 


- 225.0 

- 200 .0 


75 .0 - 125 .0 - 75 -0 - 25.0 

- 150.0 - 100 ■{) - 50*0 

REAL PRRI 


0 -0 


Figure 82. Root loci for eigenvalues No. 3 and No. 4 of the nominal HPFTP linear model. 
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Figure 83. Imaginary part versus speed for eigenvalues No. 3 and No. 4 of the 

nominal HPFTP linear model. 
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Figure 84. Critical damping ratio versus speed for eigenvalues No. 3 and No. 4 of the 

nominal HPFTP linear model. 
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Figure 85. Eigenvector display for eigenvalue No. 3 at 40,000 r/min for the 

nominal HPFTP linear model. 





While systems of this sort are generally less stable at higher speeds, the interaction 
between eigenvalue loci creates a range of speed where stability is actually enhanced. This 
particular occurrence appears to be related to the asymmetry in the housing, in particular, 
asymmetry in the frequency-dependent dynamic impedance of the housing which is an integral 
part of the rotor support. The stabilizing capacity of asymmetry is well known and documented 
(see references 12 and 18, for example). The housing impedance asymmetry effect can be shown 
by modifying the housing model so that the rotor dynamics do not "tune" with the housing 
dynamics and create this behavior. Figure 86 shows the root loci for the case where the 
frequencies of the third and fourth housing input modes are increased 20 percent. For this case, 
the switching does not occur, and the locus for eigenvalue No. 4 progresses fairly quickly into the 
right half plane. The stabilizing capacity of the switching behavior will complicate comparisons ol 
the HPFTP model's behavior with the simplified model results. However, the model will not be 
adjusted to remove this characteristic since it may truly be representative of the behavior of the 
HPFTP. 



REAL PARI 

Figure 86. Root loci for eigenvalues No. 3 and No. 4 of the nominal HPFTP linear model 

with modified housing modal data. 

The linear analysis results for the nominal HPFTP model indicate that the model is too 
stable to exhibit limit cycle instability. Based on the analytical results from the simplified model 
there are two ways to reduce the stability of the system: increase the cross-coupled stiffness 
parameter a and reduce the pertinent natural frequency (most likely by reducing the rotor 
support stiffness). The first method will increase the frequency ratio of a limit cycle while the 
second will not. The specific test data being examined here indicates a limit cycle frequency ratio 
of 47 percent; therefore, it is undesirable to adjust the model by increasing the parameter a. The 
stiffness of the bearings and seals was decreased by trial and error until suitable results were 
achieved. Greater reduction was made in the seal stiffness than in the bearing stiffness in order 
to maintain a strong effect from the nonlinearity (y parameter). The final values chosen were 
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92 percent of the nominal bearing stiffness and 60 percent of the nominal seal stiffness. The 
cross-coupled stiffness parameter for the pump interstage seals (cr) was reduced from 0.6 to 
0.55 and the turbine interstage seal and Alford effect cross-coupling was reduced to 75 percent of 
the nominal value. These changes were based on results of simulation trials. 

The linear analysis results for the modified model are presented in figures 87 through 89. 
Figures 87 (root loci) and 89 (critical damping ratio) show the same switching behavior as the 
nominal model. The transition occurs at a higher speed for the modified model. The frequency of 
eigenvalue No. 3 is about 310 Hz at the 109-percent power level speed. The onset speed of 
instability is 44,000 r/min. At the imaginary axis crossing, the frequency of the eigenvalue is 
approximately 328 Hz. This translates into a ratio of 45 percent. These values compare favorably 
with the guideline values inferred from the test data. 



Figure 87. Root loci for eigenvalues No. 3 and No. 4 of the modified HPFTP linear model. 
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D. Nonlinear Simulation Results 


The general behavior of the simplified model was characterized in terms of its behavior 
under certain restricted conditions in section V. The restricted conditions were analyzed in 
sections III and IV. The homogeneous equilibrium, unbalance mass equilibrium, and side-force 
equilibrium and their respective linearizations provided bounds for parameter ranges where limit 
cycle behavior is possible. The subharmonic response harmonic balance solution provided bounds 
for parameter ranges where subharmonic response is possible. These results were based on 
analytical and combined analytical/numerical solutions to the nonlinear system equations. These 
approaches were rather straightforward for the simplified model. Analogous methods for the more 
general and complex turbopump model have not been developed; however, the approach of using 
the restricted case results to characterize the general case results is still valid. Simulation must 
be used rather than analytical means to obtain the restricted case results. The HPFTP model 
will, therefore, be simulated for the three restricted excitation cases analyzed for the simplified 
model. The first is the homogeneous case. This case is not truly homogeneous since noise 
excitation is imposed, but it will serve to characterize the lower limit, amplitude, frequency, and 
global stability limit of the limit cycle behavior. The second is the unbalance mass equilibrium 
case (again with noise excitation). This case will show the limits where stable unbalance 
equilibrium is possible and where limit cycle will exist along with the unbalance response. The 
third is the side-force equilibrium case (also with noise). This case will show the limits where 
limit cycle is possible under the stabilizing influence of side force. These cases will be examined 
by executing a speed ramp up to just beyond the global onset speed of instability and then back 
down. A general loading case will then be examined with the same ramp profile. The response to 
this loading can then be characterized in terms of the limits and characteristics of the three 
restricted cases. The subharmonic response analysis of the simplified model was not based on a 
restricted excitation but, rather, on a specific assumed form of solution. No special simulation 
case is required to investigate this behavior. The general excitation case will be examined to 
determine whether subharmonic response occurs. 

Simulation results for the homogeneous case are presented in figures 90 and 91. Figure 90 
displays the V axis relative displacement at the inboard turbine end bearing location. The top 
graph shows the response to the ramp up, and the bottom graph shows the response to the ramp 
down. The speed profile consisted of a 4.5-s ramp from zero to 45,000 r/min and a 2.5-s ramp 
back down to 20,000 r/min. Figure 91 displays the cascade spectral plots corresponding to figure 
90. Limit cycle instability initiates at approximately 25,000 r/min with a low amplitude. The 
frequency ratio is about 55 percent. The amplitude increases steadily until a speed of about 
32,000 r/min is reached. Between this speed and about 35,000 r/min, the amplitude decreases 
slightly and the frequency ratio drops to about 47 percent. This transition corresponds to the 
eigenvalue switching observed in the linear analysis. As speed is increased beyond this 
transition, the amplitude increases in the same manner as the limit cycle in the simplified model. 
The reverse behavior occurs on the down ramp with no apparent hysteresis. The mass unbalance 
excitation case produced similar results (figs. 92 and 93). The limit cycle initiation was 
suppressed until about 32,000 r/min by the unbalance equilibrium. On the down ramp, the limit 
cycle was maintained until about 29,000 r/min. This is similar to the hysteretic behavior observed 
for the simplified model. 

The side-force excitation exhibited the greatest amount of hysteresis (fig. 94). On the up 
ramp, the limit cycle did not initiate until the global onset speed of instability was reached 
(44,000 r/min). On the down ramp, the limit cycle was sustained until about 35,500 r/min. 
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Figure 90. Homogeneous simulation results for modified H 
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Figure 93. Cascade spectral plots corresponding to 92. 
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Figure 94. Nominal side-force excitation simulation results for modified HPFTP 
nonlinear model. Inboard turbine bearing response. 


This value represents the lowest speed for which a limit cycle is possible with the given side- 
force excitation. This limit must be lowered in order to obtain general loading results which 
exhibit limit cycle at the 109-percent power level speed (36,300). This was accomplished by 
reducing the side force to 75 percent of its nominal value. This produced the results shown in 
figures 95 and 96. The limit cycle was sustained in this case until about 31,000 r/min. 

Results have been presented for the three restricted cases with bounds identified for the 
occurrence of limit cycle instability. The general combined loading case will now be presented. 
The initial case simulated consisted of a direct combination of the previous cases, i.e., nominal 
mass unbalance, 75 percent of the nominal side force, and the same noise excitation. This case 
exhibited limit-cycle behavior but did not exhibit any subharmonic response. One of the 
objectives of this examination is to determine if the behavior of the simplified model 
(subharmonic entrainment and transition to limit cycle) extends to the more complex HPFTP 
model. In order to determine this, the mass unbalance was varied by trial and error in an attempt 
to find a combination of parameters which resulted in subharmonic. Increasing the unbalance 
uniformly by 50 percent produced the desired results. 

Figures 97 and 98 present the time and spectral data for the response to the same profile 
as the three restricted cases. As speed reaches about 30,000 r/min, limit cycle initiates at a low 
amplitude with a frequency ratio just above 50 percent. At about 33,000 r/min, this limit cycle 
becomes entrained by subharmonic response at the 50-percent ratio. This entrained response 
transitions to a limit cycle of about 47-percent frequency above 35,000 r/min. The down ramp did 
not exhibit subharmonic entrainment but behaved very similarly to the unbalance excitation case 
(figs. 92 and 93). The up ramp behavior can be more clearly seen in figures 99 and 100. This case 


99 







RAMP UP RELATIVE Z DISPLAC 



OOOUD^rCNOOsJ^rCOCOO OODUD ^CSlOCsI^lDOOO 

— <000000000-^ — OOOOOOOOO— ' 

ooooooooooo ooaoooooooo 


ooooooooooo 

I I i I I 
- Z U I U 03 


ooooooooooo 

I I I I I 
• — 1 Z L) I Id 03 


101 


Figure 97. Combined excitation simulation results for modified 
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is identical to the other except that the speed profile was altered about 32,000 r/min. For this 
case, speed was increased from 32,000 r/min to 37,000 r/min in 4.5 s, the same length of time as 
the original ramp from zero to 45,000 r/min. The three zones of behavior are observable in both 
the time data and in the spectral data. The movement from the higher frequency ratio to the lower 
frequency ratio was observed for the homogeneous case and the mass unbalance case and was 
attributed to the eigenvalue loci switching behavior. However, the entrainment at the 
subharmonic frequency is unique to the general combined loading case. 

Another important result which was presented for the simplified model in section V (figs. 
76 and 77) is the nonuniqueness of the solution. This behavior is demonstrated for the HPFTP 
model by holding the speed profile for this case at 35,000 r/min. The time response for this 
simulation is shown in figure 101 with an expanded portion shown in figure 102. The spectral 
results are shown in figure 103. The system is clearly jumping between limit cycle entrained by 
subharmonic response and pure limit cycle. The frequency of the limit cycle at this speed is close 
to 0.5 due to the shift from the higher to lower ratio noted earlier. This proximity and the rapid 
manner in which the transitions occur make it difficult to obtain FFT results as clear as for the 
simplified model. The only change occurring in the parameters of the model is the random noise 
excitation. The response remains in one form for as long as 0.3 s. 

The results of the combined excitation case have clear implications related to 
interpretation of test results. For one set of parameter values, the model exhibits subharmonic 
response, limit cycle with a frequency greater than subharmonic, and limit cycle with frequency 
less than subharmonic, all within a narrow speed range. The model also exhibited transitions 
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between subharmonic and limit cycle at a fixed speed. All of these behaviors have been observed 
in test data. The model results indicate that the specific ranges of occurrence, frequency ratios, 
amplitudes, and transition speeds are quite sensitive to many model parameters which are not 
known with much certainty. The nonuniqueness and sensitivity of the solutions can easily 
account for the behaviors observed in test data. 

The results of the previous sections for the simplified model have been shown to extend 
to a more complex and realistic model of the HPFTP of the SSME. This was accomplished by 
examining engine test data, performing linear analysis of the model, and simulating the model for 
various conditions. The conclusions of this and the previous sections will be summarized in the 
next and final section. 


VII. CONCLUSIONS 


An extensive investigation has been conducted of the interaction between limit-cycle 
instability and subharmonic response in a rotordynamic system. The primary tool used in the 
study was a dimensionless, normalized model of a single mass rotor. Equilibria were determined 
for various excitations; linearizations and stability analyses were performed. A harmonic balance 
procedure was implemented to analyze subharmonic response potential. The model was 
simulated for the conditions which were analyzed and then for conditions which could not be 
treated analytical. Generalization of the results to a complex, realistic model were confirmed by 
examining the HPFTP of the SSME. This was accomplished using linear analysis and nonlinear 
simulation and by examining engine test data. The analyses and simulations were conducted 
using a general turbopump rotordynamic analysis package developed for this research. The 
conclusions which are drawn from this work can be separated into two groups: conclusions 
regarding the characteristics of the behavior, and conclusions regarding analysis and simulation 
methods. The conclusions are summarized below. 


A. Characteristics of the Behavior 

The most significant conclusion from this work is the determination that subharmonic 
response can entrain self-excited limit-cycle oscillations in rotordynamic systems. There is an 
important implication from this conclusion. The occurrence of a subsynchronous vibration at a 
frequency exactly equal to one half the shaft rotational speed is sometimes interpreted as 
evidence that the vibration is only subharmonic response and is benign in the sense of instability 
(subharmonic response may have adverse results due to overloading and fatigue which are 
unrelated to stability). These results indicate that this interpretation should not be made in any 
system in which self-excited vibration is possible. This covers virtually all turbomachinery which 
operates above a system critical speed characterized primarily by rotor motion. 

Another conclusion drawn from this investigation is that the behavior under given 
conditions is nonunique. This is known for nonlinear systems in general; however, the specific 
possible solutions were determined and demonstrated. The most striking demonstration was the 
repeated transitions between limit-cycle instability and subharmonic entrainment in the presence 
of random noise excitation. 
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Finally, the ability of side-force and mass-unbalance excitations to suppress the limit- 
cycle behavior was demonstrated. The side-force excitation may, in addition, stabilize the 
system in the small for speeds beyond the global onset speed of instability. 


B. Analysis and Simulation Methods 

Nonlinear simulation studies typically consist of extensive studies involving numerous 
variations in model parameters and excitation forces. These studies begin with the best available 
estimates of these parameters and forces. One beneficial approach identified in this study is to 
simulate the model with no excitation (homogeneous case) except random noise. The results of 
this simulation provide a bound on the speed range in which limit cycle is possible with the 
excitations applied. In addition, the amplitude for this case bounds the limit-cycle amplitudes, 
and the frequency is near the frequency obtained under general excitation. These bounds provide 
a basis for determining whether extensive investigation of limit cycle is needed. 

The important benefit of applying representative random noise excitation was clearly 
evident in this study. Without perturbation, many significant occurrences of limit cycle will be 
missed. It is also important for the perturbation to be realistic so that conclusions regarding 
comparisons with test data will not be erroneous. 

It became evident during the conduct of this study that simulations must be conducted for 
representative time durations. The transient growth of a limit cycle or the transition between 
subharmonic entrainment and limit cycle might easily be missed with short durations. 


C. Suggested Future Research 

While many important conclusions were reached in this investigation, there remain many 
unanswered questions and opportunities for further work. As with any analytical and numerical 
investigation, it is important to obtain experimental verification. The engine test data provide 
some verification for this work; however, a laboratory experimental program with well-defined 
conditions and adequate instrumentation would provide much better verification of the 
fundamental assumptions and conclusions. With regard to the HPFTP model, improvements in 
the certainty of the model parameters would simplify investigations and correlation with test 
data. Even with better parameter data, however, absolute certainty in parameter values is 
impossible. Therefore, an extensive parametric study of the HPFTP conducted in light of the 
results of this research would provide important new insight into the characteristics of the actual 
machine. This should include a penetrating review of the extensive data base of engine tests. 

With regard to test data interpretation, although the threat of divergent instability cannot 
be ruled out when subharmonic entrainment occurs, it may be possible to deduce a qualitative 
margin of stability from its occurrence. The subharmonic response analyses conducted for the 
simplified model indicated that the subharmonic ceased to exist above a threshold speed that 
was near, but below, twice the linear (zero deadband assumption) system resonance. The 
specific value depended on the system parameters. If subharmonic occurs in the response, the 
system must be operating below twice this resonance. If the destabilizing forces are known to be 
characterized by frequency ratios less than one half, then the global onset speed is known to be 
greater than twice this resonance. An investigation should be made into the potential for 
developing such an indicator of stability. Development of analytical tools for the complex model 
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which are analogous to the homogeneous equilibrium and subharmonic response harmonic 
balance methods applied to the simplified model would facilitate an investigation of this. These 
tools would provide characteristics and ranges of occurrence of the limit cycle and the 
subharmonic response without requiring costly simulations. 

Finally, a word about the basic modeling assumptions made in this work. The assumption 
was made here that the rotor was supported by a combination of linear and nonlinear support 
elements. The nonlinear support elements were treated as piecewise linear elements in radial 
deflection with only two linear regions (i.e., one break point in the load deflection curve). One 
significant situation that cannot be treated using this assumption is that of rotor-stator rubbing. 
A model of rubbing interaction by itself is very similar to a bearing with deadband. The stiffness 
and the clearance are usually larger than for a bearing and an addition force is included; the 
tangential friction force. This interaction could be handled with the bilinear assumption with the 
appropriate additions for the friction. However, the difficulty arises when there are bearings with 
deadband and rub with a much larger clearance (or deadband). This combination requires 
multilinear approximations. In a realistic model such as the HPFTP, these effects would be 
distributed as well. The implications of this for the results and conclusions presented here would 
be to introduce an additional level of potential equilibria and limit cycle instability. For example, 
for speeds beyond the global onset speed of instability in the current model, the amplitude of the 
vibration would grow until the rotor began to rub on the housing. This would, then, create a new 
limit cycle at a higher amplitude with its own global stability limit. While general results can be 
inferred by scaling the results already obtained, the application of these results to a realistic 
model such as the HPFTP require careful scrutiny. The HPFTP model examined here did not 
include the rub model and it produced results similar to test data. However, this might also be 
accomplished using a very different "recipe" of parameters (representing a less stable condition) 
which includes the rub model and its amplitude limiting effects. In order to assure that erroneous 
conclusions are not drawn regarding the HPFTP, or any machinery that is known to rub, the 
parametric study discussed above should include models of rubbing interaction. This should be 
anchored to reality by closely examining the wear which occurs in the test hardware due to 
rubbing. The mechanical work done by the rub which occurs in the model should be compared with 
the observed wear. This requires the development of material wear models which relate the work 
to estimated wear. In addition, the power loss due to the rub should be compared to total 
turbopump power and its impact assessed using an overall engine system model. This would 
determine whether the power loss would be observed in the engine performance and, if so, the 
test data related to engine performance should be examined. These reality checks (wear and 
power loss) would provide the guidance needed to determine which assumptions and parameter 
"recipe" most accurately represent the actual machinery. 
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APPENDIX A 

TURBOPUMP MODEL DEVELOPMENT 


Introduction 

A general turbopump model has been developed in order to carry out the numerical 
studies in this investigation. This model has been implemented in a package of FORTRAN 
computer programs that is referred to as the Turbomachinery Rotordynamics Analysis 
Package. There are three basic components to the package: a linear eigenvalue analysis 
provides stability and critical speed information; a linear forced response analysis provides 
steady-state response to static and dynamic loads; and a nonlinear time domain simula- 
tion provides the total solution (transient and steady-state) and incorporates important 
non-linear effects such as bearing clearance (deadband) and seal rubbing. The simpli- 
fied model developed in chapter II can be implemented using this package by specifying 
the input data appropriately. This appendix covers the development of the equations of 
motion for the model and the solution procedures employed. 

Conceptual Model 

The turbomachinery rotordynamics analysis package is based on a conceptual model 
of a symmetric flexible rotor supported in a nonsymmetric flexible housing by flexible 
connection elements. This is illustrated schematically in figure 104. The rotor is charac- 
terized by its free-free normal modes of vibration. Likewise, the housing is characterized 
by its free interface (no rotor) normal modes. One axial rigid body degree of freedom 
is included for the rotor. This is included to couple with the axial component that may 
be present for each housing mode. This coupling usually takes place across a hydrody- 
namic thrust balance piston. Damping is added to the rotor and housing by specifying a 
damping ratio for each mode. Gyroscopic effects are included as generalized forces on the 
rotor and they create coupling between the rotor modes. The normal modes of the rotor 
and housing must be predetermined using structural models or other available means. 
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Figure 104. Turbopump conceptual model schematic. 


The connection elements are typically rolling element bearings and fluidic seals and are 
characterized by stiffness, damping, and inertia matrices that can contain coupling terms 
between orthogonal lateral axes (cross-coupling). In general, the elements of these ma- 
trices are functions of engine power level. This is caused by such factors as centrifugal 
loading on ball bearings and different pressure drops across fluidic seals. Since turbopump 
rotational speed is also a function of power level the stiffness, damping and inertia co- 
efficients can be expressed as functions of speed. These functions can be represented in 
a number of ways. In its present form the package can use two methods to define each 
function: a polynominal in speed, or a table lookup. The polynomial coefficients must 
be pre-determined by curve fitting tables of data for each function. The table lookup can 
use linear or Hermite cubic interpolation. 

The equations of motion for the turbopump system are developed by deriving the 
equations of motion for the rotor and housing separately. The forces due to the relative 
motion across the connection elements are then added as generalized forces acting on 
the rotor and housing. The coordinate system used to define the model is an inertial, 
right handed system with the x axis along the undeflected rotor centerline. The y and z 
axes are in orthogonal, lateral directions. The orientation of the y and z axes is usually 
determined by the structural model of the housing since the rotor is symmetric. Care 
must be taken to ensure that the proper algebraic sign is used for the rotation speed 
based on the right hand rule for the coordinate system. 

Rotor Equations of Motion 

The rotor is treated as a collection of rigid bodies. The equations of motion for 
the rotor are derived using Lagrange’s equations. The kinetic energy, potential function, 
dissipation function, and virtual work expression are first defined in terms of the physical 
coordinates of the individual rigid bodies. These functions are then expressed in terms 
of the rotor free-free normal modal coordinates. The equations of motion are obtained 
by substituting these functions into Lagrange’s equations. 

The coordinate system used to define the motion of an individual rigid body segment 
of the rotor is shown in figure 105. The x, y, z axes define the inertial reference frame. 
The x, y, z axes define the body fixed reference frame. The Euler angles 9 y , 9 Z , and 9 X 
are defined in the figure. 
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The kinetic energy for the i th rigid body can be written as the sum of the translational 
and rotational kinetic energies 


T i = im i (x? + j? + i?)+5 
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where u> £ , u> 5 , and u ; 2 are the body axis components of the angular velocity. Since the 
x axis is taken to be the axis of symmetry and each body is considered to be a body of 
revolution, the product of inertia terms will be zero and equation (Al) can be written as 


Ti = + y} + if) + \la + ^(wj. + w l ,) 


(A2) 


The angular velocity components u > tt , a> 5 , , and u- z< can be written in terms of the Euler 
angle rates as follows 


^x, = Ox, + 6y. sin (0 Zi ) 

u>y t = 9 yi cos (6 Zt ) cos (6 Xi ) + 0 Zx sin (6 Xt ) (A3) 

u Si = -8 Vi cos(6 Zi ) sin(0x, ) + 0 S . cos (0 Xi ) 

Substituting equation (A3) into equation (A2) yields 


Ti = \mi{x] + y \ f + z\ ?) + \l a , [9 X , + 0 Vi sin(0 r ,)) 2 + \h,{0 2 yt cos 2 (6 Zi ) + 0%) (A4) 


The rotor will be constrained by the bearings to have small motion, therefore, the angles 
Q y . and 6 Zi can be considered to be small quantities. Substituting Taylor series expansions 
for sin(0 2i ) and cos(# 2( ) into equation (A4) yields 


Ti = 
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Since the kinetic energy expression will be differentiated once in Lagrange’s equations, 
terms of second order in small quantities must be retained here. Neglecting terms of 
higher order yields 


T, = + y? + z}) + + *1 ) (A6) 

The total kinetic energy of the rotor is the sum of the individual kinetic energies, 


r = 2> (at) 

i 

It is assumed here that the rotor is axially and torsionally rigid so that equation (A7) 
becomes 


T 


+ 5E / ^+»«> +n (E / ‘AA-) 

t ' » ' 


where Q replaces 6 Xi . Equation (A8) can be written in matrix notation as 


(A8) 


where 


T = 


-Mi 2 + -n T I a n + ^y T my + |z T mz 
+ l -Q?i t Qy + \eji t @ z + a©^i & e z 


(A9) 


M = ^Tm, 

i 


(A10) 


y = 


' 2/1 

V2 

< 

“ Vn 


(All) 
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and similarly for z, © y , and © z , 



‘ mi 0 ... 0 

0 mo 0 

. 0 0 ... m n _ 


(A12) 


(A13) 


and similarly for I t and I a . 

The potential and dissipation functions will be defined in terms of the rotor physical 
coordinates. Since the rotor is axially and torsionally rigid, only the y, z, © y , and © 2 
coordinates will contribute to these functions. The potential function is written 



k w = 


1 

2 


„ T 

y 1 

K yy 

Kyt. 

0 

0 - 

( y 'i 

©z 1 

K 6*y 

K B,6, 

0 

0 

1 ©z 

z 

0 

0 

Hzz 


1 2 1 

©y ) 

. 0 

0 

K6 v z 

- 

l© y J 


(A14) 


where k represents the rotor stiffness matrix; the upper left and lower right quadrants 
are identical (symmetric rotor). 

The dissipation function can be more easily treated using coordinates that rotate 
with the rotor at speed ST. The following transformation defines the new coordinates 


r < 


cos(Qf) sin(Qf) 

/ Wy \ 

— sin(flt) cos (fit) 

IwJ 


(A15) 


where 


w * = { e, } (A16) 

w * = {e,} (A17) 
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(A18) 


In terms of the rotating coordinates r, the dissipation function is written as 


G = -r T B r 
2 

Differentiating equation (A15) with respect to time gives 


_ [ cos (fit) sin(Dt)l f w y + ftw 


[ — sin(Df) cos(flt)J \ w z — fiw y 
Substituting equation (A21) into equation (A20) yields 


Q _ M Wy + flW z 
U 2 I W Z - flWy 


cos(fii) — sin(fit) 
sin(Df) cos(Df) 


B v 0 
0 B ( 


cos (fit) sin(fif) f w y -f flw z 
- sin(flt) cos(f It) \ w z — flwy 


(A22) 


where B,, and B^ are identical and represent the damping matrix for the nonspinning 
rotor. Equation (A22) can be rewritten as 


r y + flz 

1 G) Z + flQy 

G= 2 z - fty 

k © y - n© z 


y + flz 

© Z + fl©y 

z — fly 
0y - n© z 


(A23) 


The generalized forces acting on the rotor are treated through the virtual work 
expression. The virtual work is the product of a virtual displacement of a coordinate and 
the component of a generalized force acting on the coordinate. The virtual work can be 


written as 


&W = Y^[( F x,^x, + F Vi e yi + F Zi e Zi ) • (6xii x> + %e v , + 6zji Zt ) 


+ (T Xi e Xt + T Vt iy t + T Zt e Zi ) ■ (SO x .e Xi + 69 y> e yi + 86 Zx e z ^)\ (A24) 
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where e q is a unit vector along the q axis. Equation (A24) can be simplified to 


SW = + F Vi 6yi + F Zi 6zi + T Xi 88 Xi + T yi 68 Vi + T Xi 88 yi (e £i ■ e Vi ) 

i 

+ Ty t 68 Xi (e yi • e$ ( . ) + T Zi 88 Xi (e Zi • c*. ) + T Zi 68 Zi (e z< - e 2 ')] (A25) 


From figure 105, the unit vector dot products are determined to be 


e Xi • e y , - e Vl • in = sin(0 2i ) « 8 Zi ( A26) 

e Zi •«*, = - sin(0 y , ) « -0 y , (A27) 

e 2i • e 2 < = cos(# yi ) « 1 (A28) 

for small angles. Substituting equations (A26) - (A28) into (A25) and simplifying yields 


6W = £[ F m 6xi + F yi 6yi + F Zi 6zi + (T Xi - T Zi 0 yi + T y ,8 Zt )68 x> 

i 

+ (T„e, i +T, l )M, l +T,.M K } (A29) 

The torques T Vx and T Zi will be linear functions of the displacements 8 y% and t 9 Zt . Since 
these angles are considered to be small quantities, the products of the torques T Vi and 
T Zi and these angles become proportional to terms of second order in small quantities. 
Neglecting terms of second order and higher in small quantities simplifies equation (A29) 
to 


sw = + F Vi 6yi + F Zt Szi + T X M X , + T z ,88 Zt + (T y . + T Xt 8 z ,)68 y ,] (A30) 


The kinetic energy, potential function, dissipation function, and virtual work have 
all been defined in terms of the rotor physical coordinates y, z, © y , and © z . These 
functions can be expressed in terms of the rotor free-free modal coordinates using the 
following transformations 


y = $q y 

(A31) 

© z = ^q y 

(A32) 

z = $q z 

(A33) 

© y = -*q z 

(A34) 


where # is the rotor free-free modal displacement matrix, 'I' is the rotor free-free modal 
rotation matrix, q y is the rotor free-free modal coordinate vector in the y — x plane, and 
q z is the rotor free-free modal coordinate vector in the z - x plane. It can be noted that 
the rigid body modes must be included in $ and Substituting equations (A31)-(A34) 
into equation (A9) yields the kinetic energy in terms of rotor modal coordinates: 


T = 


The rotor free-free 


+ ^ft T I a ft -I- ^qJ# T m#q y + iqJ$ T m$q z 
+ iqJ* T I t *q y + ] -qJ* T l t *q x - ClqJ* T I & *q y (A35) 

normal modal vectors are orthogonal and properly scaled so that 


$ T m$ + $' T I t \I' = I 

Therefore, equation (A35) reduces to 


(A36) 


T = ~Mx 2 + ^n T I a ft + ^q y q y + ^qjq z - ftq^Tq 


(A37) 


where 


T = ^ T I a ^ (A38) 

Substituting equations (A31)-(A34) into equation (A14) yields the potential function in 
terms of rotor modal coordinates: 
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[ $qy 'j T f K yy K yB, 0 0 1 ( 

T/ - I i ^ qy 1 Ke ‘ y Ke ‘ e ‘ 0 0 ) *q y 

2 | $q z 1.0 0 K zz K z e y | #q* 

( -tyq z ) 0 0 K 6yZ *e y 6 y - I -^q z 

Due to the orthogonality and proper scaling 


r*T ^Ti K yy K y8* ® =u 2 (A40) 

1 l K e, y Ke.e, J L'*' J 

and 



which gives 





where 



>, t r 2 

q y 1 w n 

q z J 0 




(A41) 


(A42) 



(A43) 


Substituting equations (A31)-(A34) into equation (A23) yields the dissipation function 
in terms of rotor modal coordinates: 


' #(q y + Dq z ) ^ T Bt,^ 0 0 j | $(q y + Dq z ) 'j 

lj ^(q y -Dq z ) I Be (X , Be ( e ( 0 0 I ^(q y -Dq z ) I 

^ — 2 | #(q z - fiq y ) ( 0 o B cc B ce>| I $(q 2 ~ Oq y ) 

, -*(q z + Dq y ) J L 0 0 B*„ ( B m ,J l -^(q* + ^y) > 

(A44) 

If the rotor damping is assumed to be proportional damping, the damping matrix can be 
diagonalized by the modal matrices just as the stiffness matrix was diagonalized. This 
allows equation (A44) to be written as 
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2 C w n 


2Ci u \ 


Ui 0 

0 2C 2 w 


n a 


0 - 
0 


2Gtw nt 


(A46) 


Substituting equations (A31)-(A34) into equation (A30) yields the virtual work expression 
in terms of rotor modal coordinates: 


6W = F x 6x+6QjT x + 6q?($ T F y + * r T z )+6q? [* T F Z - * T (T y + T^q y )] (A47) 
where 


T* - 

x X 



L 0 


0 ... 0 - 
T X2 ... 0 

0 ... T Xk . 


(A48) 


is the diagonal matrix formed from the vector T x . 

The kinetic energy, potential function, dissipation function, and virtual work expres- 
sion are now in the form desired for use with Lagrange’s equations. Lagrange’s equations 
can be written in the following form 


dt V dq { ) dq t + dqi + dq { 


(A49) 


where is the generalized force acting at coordinate <7, and is implicitly defined by 


6W = Y / Qi^i 

t 


(A50) 
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Substituting equations (A37), (A42), and (A45) into equation (A49) and comparing equa- 
tion (A47) with equation (A50) yields the following set of equations of motion for the 
rotor: 


Mx = F x (A51) 

I a H - qjTq y - q?Tq y = T x (A52) 

q y + £lTq z + u^q y + 2<> n q y + fi2f>„q z = * T F y + * T T z (A53) 

q z -Q,rq y + ulq z + 2(;u n q z -n2(u n qy-nTqy = $ T F Z - * T T y - ^ T T^q y (A54) 
Neglecting second order terms for small angles, it can be seen from equation (A52) that 


T x = I a 17 (A55) 

However, since the rotor is torsionally rigid, all elements of the vector 17 are equal, and 
therefore 


t; = ni a (A 56 ) 

Substituting equation (A56) into equation (A54) yields 

q z - firq y + u 2 n q z + 2(w n q z - 172Cw„q y = * T F Z - * T T y (A57) 

Rewriting the complete set 

(A58) 
(A59) 
(A60) 


1*17 = T> 


x = 


M 


q y + 2(cj n q y + 17rq z + + S72£o> n qz = ^ T F y + ^ T T z 
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q z + 2(u n q z - ftTq y + u>lq z - ft2Cu> n q y = $ T F Z - ^ T T y (A61) 

Equations (A58)-(A61) are the rotor equations of motion. Equation (A58) defines 
the torsional motion. The rotor speed ft will be a specified function of time, therefore, 
equation (A58) will not be needed. The forces and torques F x , F y , F z , T y , and T z contain 
the linear and nonlinear forces due to the elements connecting the rotor to the housing 
as well as externally applied forces. 

Housing Equations of Motion 

The equations of motion for the housing are derived in the same manner as for the 
rotor with the following exceptions. The housing doesn’t rotate so that all terms involving 
ft are absent from the housing equations. Also, the modal matrix used to diagonalize the 
housing equations is derived from fixed-free boundary conditions. However, the housing- 
rotor interface coordinates are free. The housing is not symmetric and its modes are not 
entirely planar so the modal coordinates can not be separated into y and z coordinates as 
in the case with the rotor. With these exceptions noted, the housing equations of motion 
can be written directly 


P + 2C h w„ h p +wjj h p = 


Fk. + *£F hy + *£ F h , + tfjT hjr + *£T h , (A62) 


Combined System Equations of Motion 

With the rotor and housing equations determined, the system equations of motion 
can be written, 


[ * ' 

1 q y 
v .. 

► + 

q ‘ 

l P J 



0 0 0 
0 2C r u> nr ftT 

0 -ftr 2 £ r w„ r 

0 0 0 


0 ■ 


' X ' 


'0 

0 

0 

0 • 


X ) 

0 

< 

q y 

► + 

0 

CJ 2 

n r 

ft 2 CrW Dr 

0 

i 

qy l 

1 

0 


q z 

0 

— n 2 ^ r u? nr 

i J 2 

n r 

0 


q z 

^Ch^ni, - 


, p > 


.0 

0 

0 

u >2 

n h J 


l p J 


/ 


M 


$ T F y + * t T z 
* t F z - ¥ T T y 




hi 


+ *£Fh, 


+ *£f„. + 


*£Th y + * h T ,T h , 


(A63) 
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At this point the rotor and housing appear to be uncoupled. However, the force vectors 
on the right hand side of equation (A63) are partially due to rotor-housing interaction 
and can be expressed as functions of the modal coordinates. First, the forces are written 
in terms of physical coordinates: 

ijr -E/i ) C x (i i/i) 

(A64) 

h 

II 

1 

(A65) 

where k x and c x are the axial stiffness and damping coefficients connecting the rotor to 
the housing; 

F y = -k y y - Q y z + k y yi, + Q y z h - c y y - C Qy z + c y y h + C Qy z h 
- m y y - M Qy z + m y yi, + MQ y zu + FEr y + F„ y 

(A66) 

F z = Q z y - k 2 z - Q z yh + k z z h + CQ,y - c z z - C Qi yh + c z z h 
+ Mqj - m z z - MQ,yh + m z z h + Fcr, + F n , 

(A67) 

Fh y = — (F y — F Ery ) ~ F E hy 

(A68) 

Fu, — “(F z — F Ers ) — F E h s 

(A69) 

where the coefficient matrices k, Q, c, Cq, m, and Mq, are diagonal; 


T z = — kt © z 4* k t 0 Zh 

(A70) 

Ty = — k t ©y + kt©y h 

(A71) 

T hl = — T z 

(A72) 

Th y = -T y 

(A73) 
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where the coefficient matrix k t is diagonal. The external excitation forces (F Er and Feu) 
and nonlinear interaction forces (F„) will be defined in later sections. 

The physical coordinates (x, y, z , etc.) can be replaced by the modal transformation 
given by equations (A31)-(A34) for the rotor and the following for the housing: 


*h = $h x P 

(A74) 

yh = ^h y p 

(A75) 

©Z h = *h,P 

(A76) 

Zh = #h,P 

(A77) 

©y>. = ^h y P 

(A78) 


Performing these transformations and substituting the resulting force expressions into 
equation (A63) yields: 


X = -^(x - 4» hx p) - ^(x - # hx p) (A79) 

cjy -f- 2£ r L> nr q y + ftTq z + w„ r q y T fl2£ r u> nr q z 

= -# T ky$q y - $ T Qy$q z + ^ T ky# hy p + S^QyS^P 

- $ T c y $q y - 4> T C Qy $q z + # T Cy# hy p + $ T C Qy $ hl P 

- $ T m y $q y - $ T M Qy <i>q z + <& T m y $ hy p + # T M Qy # hi p 

- $' T k t 'I'q y + ^ T k t ^ h ,p + # T (F Ery + F ny ) (A80) 

T 2£ r u> nr q z — nrq y + w ny q z — f!2£ r tc> nr q y 

= $ T Q z #q y - $ T k z 4>q z - 4> T Q z $ hy p + $ T k z # hi p 
+ $ T C Qj $q y - $ T c z $q z - * T C Q ,* hy p + # T c z $ hi p 
+ ^MQ^qy - $ T m z #q z - S^MQ^UyP + 3> T m z 3» hl P 

- ^ T k t ^q z - ^ T k t ^ hy p + * T (F Er , + F n J (A81) 
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p + 2CbU>n h P + Wn h P = ~ + &b x C * X ~ ^h x C ^h x P 

+ *i7 y k y*qy + *jQy*q* ~ *£ k y*h,P - *jQy*h,P 

- #£Q z $q y + $Jk z #q z + $£Qz$h y P - $£ k z$h,P 

+ ^hy c y ^q y + *JC Qy #q, - *£c y *b y P - ^if y CQ y ^h.P 

- C Q ,$q y + *£ c **q» + *£ C Q,^h y P “ *h, c z # hl p 

+ #Jm y $q y + #£M Qy #q z - #£m y $ hy P - ^jMQ y ^h,P 

- $£M Qi $q y + $£m z $q z + $£M Qm 3>h y P - &h, m z$h,P 
+ ^ Jkt^q y - ^Jkt^h.P - ^5 k t^qz - 'J'h y k t 'I'h y P 

- * £(F B h y + F„ y ) - *£(F Eh , + F n J (A82) 


Each force expression on the right-hand side of equations (A80)-(A82) should be recog- 
nized as the sum of the generalized forces due to the physical forces, i.e., 

$ T F y = ^^F yi (A83) 

i= 1 

where <f> f represents the row of 4>. These equations can be written more compactly 
by combining the coefficients of the generalized coordinates on the left-hand side of the 
equations as follows: 



' x ' 

< X > 


r X > 

f 0 1 

M < 

q y 

q z 

- + c 1 ? y 

1 qz 

> +K< 

Qy 

Qz 

* T (F Ery + F„,) ( 

- 1 * t (Fe„ + F„,) 


‘ p > 

l p J 


, P V 

-*£ (Feu, + F„ y ) - *£ (Fek, + F„. ) J 


(A84) 


where 
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0 

0 

I + # T m y # 

* T M Qy $ 

-* t M Q i # 

I + # T m z # 

“*£“»* "I 

( -*£m q ,s 

+ »jMo,t J 

1 - 


- $ T m y $ hy 
-* T M Qy * h , 

* T M Qi * hy 
- $ T m z #h_ 


I + ^Jni y #h jr 

+ ^M Qy # hi - #£M Qi $ h> 
+ $ h T ,m I $h 1 


r w» r ] r fir 

$ T Cy$ J + * T C Qy * 

!r 1 r 20„ r 

> T Cq,^ J +$ T < 


- $ T c 


-^h T y c y ^ 

+ *£ c q,* 


-^ y c Qy # 


-6 , £l 

^h x ^ 


- $ T c y $ hy 

- * T Cq y * h , 

« T C Q ,*h, 

- ^ T c z ^, 1 


,*h y ' 
-z$h, 4 


2ChWn h +^ h T y C y # hy 
+ *h T y Cq y *h,-* h T .C Q ,* h) 
+ *£c.*h, + 


(A86) 
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hi. 

M 


0 


0 


K = 


0 


u>: 


+ # T ky^ 

+ ^ T k t ^ J 


Q,2( r u nr 


+ 4> T Qy3> 




f"n 1 

' - S12C r u>„, ' 


+ $ T k z $ 

, - $ t Q.$ , 


, + ^ T k t ^ , 


^* h * ~M 

( - S T ky$ hy 

- # T Qy^h, 

-* T k t tfh, J 

' * T Q z *h y 

- * T k z * h , 
l + \k T k t 'i'h s 



' ~ *h T y ky* ' 


' -*£Qy*' 



+ * h x ,Qz* 


-*h T ,kz* 



, - *£k t * - 





^n h + ^ky^h. 


+ *jQy*h. ~ *£Qz*h 

+ $£ k z$h, +&u x k T$h x 

{ +*£k t *h, + ^kt^h 


(A87) 

Each of these system matrices is, in general, a function of speed since the coefficient 
matrices are functions of speed. The right-hand side of equation (A84) contains only 
externally applied forces and nonlinear interaction forces. These will be discussed next. 


External Forces 

There are four external excitations that can be considered using the rotordynamic 
analyses package: rotor mass unbalance forces, static side forces, white noise forces, and 
pulse perturbation forces. The first three represent actual forces in the turbopump. The 
pulse perturbation is a tool used to study the characteristics of the nonlinear system. 

The unbalance forces are the inertia forces due to the acceleration of the eccentric 
mass of the rotor. Since they are inertial, they are applied only to the rotor. The y and 
z axis components are ±90 degrees out of phase, depending on the direction of rotation. 
The side forces are due to non-uniform circumferential pressure distributions that occur 
in the turbines and the pump discharge volutes. These forces are applied to the rotor 
and the housing with equal magnitude, but opposite direction. The white noise forces 
represent the various random excitations that occur in the turbopump. Random pressure 
fluctuations in the turbines and pumps and external acoustic noise are examples. Noise 
can be applied to the rotor and housing with equal but opposite forces, as with the side 
forces, or to the rotor and/or housing independently. The pulse perturbation is applied 


129 


at only one axial position and is applied to the rotor only. It can be defined as a square 
pulse or a short duration oscillating force with a prescribed frequency. 

The unbalance force is defined by specifying the product of mass and eccentricity 

along with the phase for each axial location. For one location, the force is written as 

F Uy , = (ma)i6 2 cos (6 + </>,) + (ma),a sin(0 + 4>i ) (A88) 

Fuz, = ( ma)i9 2 sin(0 + 4>i) - (ma),a cos(0 + fa) (A89) 

where 



Jo Jo 


(A90) 


The side force is defined by specifying the y and z axis components for each axial location 
as quadratic polynominals in pump speed 


FSy, — CyO, + £lC y i, + fi 2 C y 2, 

(A91) 

Fsz, — C z 0 , + + £l 2 C z 2 , 

(A92) 

The white noise forces (F/w yi , F/vy j( , F/v/, yi , and F/v^J are defined by scaling a uniform 
random number sequence with range (-1,1) to a desired range. The pulse perturbation 
is either of the form 

a. t- 

II 

»> 

£ 

6(t-r)-6{t-(T + T))] 

(A93) 

r A Pz 

Fp z = 

T 

6(t-T)-6{t-(T + r))] 

(A94) 

or 



Fp y = A Py cos (w p <) [$(t - r) - S(t - (T + r))l 

(A95) 
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F Pz = A Pz sin (w p t) U(t -r)-S(t-(T + r))] 


(A96) 


where 



if t > 0 
if t < 0 


The total of all external excitation forces can now be written as 


FEr v , = FSy, + fWr v , + Fu y , + f>ikFp yt 
Fet „ = Fsz , + Ftfr,, + Fuz, + {>ikFpz k 
Fsh vl = F S y, + 

= F Sl) + 


(A97) 


(A98) 

(A99) 

(A100) 

(A101) 


Nonlinear Interaction Forces 

Rotor-housing interaction forces are functions of the physical coordinates at partic- 
ular locations. The linear interaction forces can be expressed as functions of the general- 
ized coordinates using coordinate transformations. This was done to arrive at equations 
(A79)-(A82) and (A84). The nonlinear interaction forces, however, cannot be treated in 
this way. At each instant in time, the generalized coordinates must be transformed into 
physical coordinates. The nonlinear forces are then calculated as functions of the physical 
coordinates. The physical forces are then transformed into generalized forces and applied 
to the generalized coordinates. 

There are three generic types of nonlinear force elements that can be represented 
using the rotordynamics analysis package. The first is bearing clearance or deadband. 
For this type of element, no force is produced until the relative displacement between 
the rotor and housing exceeds some specified clearance. After this clearance has been 
exceeded, the force is represented as a piecewise linear spring and damper. Figure 106 
illustrates the piecewise Unear spring force versus radial deflection. Figure 107 shows the 
relative displacement and velocity vector diagram which is used to aid in writing the y 




Figure 106. Deadband nonlinear element piecewise linear restoring force 

versus radial deflection. 

and z axis components of the nonlinear force. The element can be conceptualized as 
shown in figure 108. This figure shows a separate spring damper subelement for each 
clearance. The coefficients for the second and third of these are defined to be the changes 
in the coefficients defined in figure 106. This makes it possible to write the force equation 
as if the three parts were independent when they are actually not (figure 109). Referring 
to figure 107, the actual displacement and velocity components across a given subelement 
are 


Ve, = y - Si cos 6 

(A102) 

2 e , = z — Si sin 6 

(A103) 
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Figure 108. Conceptual representation of deadband nonlinear element. 
Heavy-lined circles represent massless rings. 
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y e , = y + <5,0 sin 0 

(A104) 


i e< = i — <5^0 cos 0 

(A105) 

where 

a Z V 

e= r * 

(A106) 

and 




r = y/l f 2 + Z 2 

(A107) 

Recognizing that 

cos 0 = - 
r 

(A108) 

and 

. - 2 
sin 6 = - 
r 

(A109) 

equations (A102)-(A105) can 

be rewritten as 



<ci 

CO 

II 

<ci 
►— * 
I 

-i 

(A110) 



(Alll) 


y e , = y + — 0-2 

T 

(A112) 


z e> = z By 

T 

(A113) 


With these defined, the force components can be written directly for each range of dis- 
placement. For r < So, 
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(A114) 


Fy= 0 

F z = 0 (A115) 

For So < r < £1 , 

Fy = -kyoVeo ~ c Vo2/e 0 (A116) 

Fz = -k Zo z eo - c Zo z eo (A117) 

For 61 < r < Si , 


Fy = ~{ky , - A: yo )y ei - (c Vj - c yo )y e] - k yo y e<J - c yo y eo (A118) 

F z — ~{k Zi — k Zo )z ex — (c 2] — c Zo )z e] — k Zo z eo — c Zo z eo ( A 119) 

For $2 < r i 

Fy = — (^y 2 “ kyi)ye 2 ~ ( c y 2 ~ c y\)Ve 2 ~ (^yj — k yo )y ei 

~ ( c y, - c yo)Ve , - k yo y eo - c yo y eo (A120) 

F z = ~{k Z2 fc 2 j )z e7 — (c 2j — c Zl )z e2 (^ 2 j k Zo )z ei 

— ( c 2 j — C 2 0 )^ei — k ZQ z eo — c ZQ z eo (A121) 

Substituting equations (A110) - (A113) for each range of displacement into equations 
(A116) - (A 121 ) and simplifying yields the following: For 6 0 < T < ^i> 

Fy = ~ ~ c vo{y+ ^ z ) (A 122 ) 

F z = -k a 0 (l - ^)*-c 2 o (z- jdy) (A123) 


For < r < S 2 , 


Fy = ~k yi 


j _ (1 _ ^ y ° <^1 -Jo 


'yi 


S 1 


y - c 


yi 




-yi 


s 1 


(A124) 
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F z = -k z 
For 62 < r, 


1 ^1 /, k Zo 61 - 60 

r 1 K Si ] 


z - c 2 


S\ (. c Zo - 6 0 


( 1 -r L£ 4 L ^H (A125) 

C Zl Oj J 


— ^V2 


— c 




V2 


s 2 


F, = -k 


22 


- C 


'-Tf'-r 


22 


k yi 62 - 

s 1 _ 

ky o^l 

0 

1 

ky 2 < 5 2 


ky. 

62 n 

_ c v» h. 

- 1 h 

_ £*> 

6\ — 6q 

c yi 

62 

c y2 


k z 1 62 — 

*1 _ 

kz 0^1 

-<0.1 

k z 2 62 


k zt 

62 

1 <"> 
N 

I 

to 

-Si 

_ £f 5 L 

61 — 60 

c * 2 

62 

c * 2 

6 2 


)0Z 


)0y] 


(A 126 ) 


(A 127 ) 


Each of these force expressions is of the form 


S S 

Fy = -k y (l “ a v~)y ~ c y(v + Py~O z ) (A 128 ) 

F z = — fc-^1 - — c z (^z - fiz-dyj (A 129 ) 

The stiffness and damping coefficients in these expressions are, in general, functions of 
speed. 

The second type of nonlinear force element that can be represented is rotor-stator 
rubbing. Rubbing is very similar to bearing clearance in that there is an abrupt stiffness 
increase when contact is made. The rubbing force element contains the additional effect 
of friction which produces a force tangential to the contact surface. The rub element 
can be conceptualized as shown in figure 110 . It can be modeled using two different 
formulations. The first is the more conventional and simpler of the two. First, the radial 
force is calculated neglecting the effect of the frictional force: 


Fr = ~MM - c) = -Mr|(l “ j^j) (A 130 ) 

Then the tangential force is determined as the product of the radial force and the coeffi- 
cient of friction: 
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F e = liFr = -fxk R |r| (1 - 


(A 131 ) 




The assumption is made here that the contact surface normal vector remains collinear with 
the radial deflection vector. This assumption may be valid for low friction or light contact. 
In general, however, the contact surface normal will be collinear with the resultant of the 
radial and tangential forces. This is illustrated graphically in figure 111. A formulation 
including this effect was developed by von Pragenau [19]. A new vector S has been 
introduced to define the surface displacement. The total force is equal to the stiffness 
multiplied by the surface displacement vector 


F = -fc/i S (A132) 

This can be expressed in terms of the contact surface normal and tangential components 
as 


F = -F N e N - fiF N e T (A133) 

The angle between the force F and the surface contact normal e N is then 


7 = tan 1 n (A134) 

The displacement vectors are redrawn for clarity in figure 112. All three sides and one 
angle are known, therefore, from elementary trigonometry, we have 


from which 


sin(7r~7) sin/? 


(A135) 


sin 0 = 7-7 sin(7r - 7) = sin 7 

r r 


Noting that 


(A136) 


sin 7 = sin(tan 1 ji) = 


\fi +/* 2 


(A137) 




Figure 112. Rub model displacement vector diagram. 


and that 


yields 


a = 7r — /? — (7T — 7) = 7 — /? 


(A138) 


c* = tan 1 fi - sin ^ 7 - 7 — ~ ) 

M Mrk/TTT?/ 


■| yTT/T 2 

The force F can now be expressed in terms of radial and tangential components. 


(A139) 


but 


F = -k R S = -k R (r - ce N ) 


(A140) 
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e/v = e r cos a - eg sin a 


(A141) 


therefore 


F = -kpi [(r — c cos a)e r + c sin aeg\ 
From equations (A139) and (A136), noting that 

cos (3 — \J 1 — sin 2 /? 

it can be shown that 


(A142) 


(A143) 


cos q — 


and 


sm a = 


\fi + M 2 


\/l +M 2 


i _ VlL + c 


l - 


| r | 2 1 +/x 2 


c 2 /i 2 


|r| 1 + // 2 




(A144) 


(A145) 


|r| 2 1 + /x 2 |r| 1 + n 2 

Substituting (A144) - (A145) into (A142) yields the radial and tangential force compo- 
nents in terms of the radial displacement r 


F = - 



(A146) 


The third type of nonlinear force element is the floating ring seal. This is a very 
complex element since it involves an additional mass suspended from the rotor and in- 
teracting with the housing through friction. Since this element was not used in any part 
of this study, the details of its formulation will not be included here. 

Eccentricities are geometric offsets from the normal centerline of the housing (stator) 
or rotor. The are not inherently nonlinear effects, but require the same transformations 
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(from generalized coordinates to physical coordinates and back) as the nonlinear force 
calculations; therefore, they are calculated along with the nonlinear forces. Their most 
pronounced effects occur when nonlinearities are present by altering the deflections re- 
quired to reach a certain threshold in the nonlinear function (i.e., exceed a seal clearance 
and rub). Their effects are included by adding the offsets to the relative deflection, 
velocity, and acceleration expressions as follows: 


Vr el = Vr - Vh + (r COS (0 + </> r ) - (A147) 

Z rel ~ 4" C r sin(0 4- 4> t ) - € ht (A148) 

Vrei = Vr ~ Vh - 0( r sin(0 + (f) r ) (A149) 

Zret = Zfi 0€ T cos(0 ^r) (A150) 

Vrei = Vt - Vh - Otr sin(^ + (j ) r ) - 0 2 ( r COS (6 + 4> r ) (A151 ) 

Zrei = z r - Zh + 0( r cos (6 + (j> r ) - ^ 2 € r sin(^ + <p r ) (A152) 


where 6 is the angular position of the rotor and for constant speed cases 

6 = (A153) 


Solution Procedure 

The equations of motion for the system have been defined in terms of generalized 
coordinates associated with the component modes of the rotor and housing. The develop- 
ment assumed that a complete set of component modes would be used to transform from 
the physical coordinates to the generalized coordinates. If a reduced or truncated set of 
modes is used (thereby reducing the order of the model), the coordinate transformation 
(equations (A31) - (A34) for the rotor, and equations (A74) - (A78) for the housing be- 
come approximations and error is introduced into the model. If a sufficient number of 
modes is retained, the magnitude of the error can be made small. The error can be further 
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reduced for a given number of retained modes by a method known as condensation. This 
method partially utilizes the modes to be eliminated rather then simply truncating them. 
The method is developed for this application in reference 20. In order to avoid possible 
influences on the results, condensation was not used in the numerical investigations of 
this study. 

The equations of motion given by equations (A84) - (A87) are a set of coupled, 
nonlinear second order differential equations. The nonlinearities in these equations are all 
contained in the rotor-housing interaction forces. These equations can be linearized in two 
ways. First, the nonlinear force expressions can be linearized about some operating state 
as described in chapter III. The state is determined by running the nonlinear transient 
simulation until a specified time is reached. Second, the nonlinearities can be neglected. 
For the bearing clearance nonlinearity, this means that the clearance is assumed to be 
either zero or infinite. The stiffness is then linear and can be handled like the other 
linear stiffnesses. For the rubbing nonlinearity, the clearance is assumed to be infinite so 
that rubbing never occurs. The linear set of equations can now be solved using standard 
techniques. 

The stability analysis program obtains the homogeneous solution to the equations of 
motion (A84): 


Mt) + C 17 + Kt? = 0 


(A154) 


where 



This can be put into first order form by defining 

'-{!} 


(A155) 


(A156) 


and writing 


0 = 0 


(A157) 


Letting 


0 

M 



0 

K 


and rearranging yields 



I 

C 



0 

K 


(A158) 

(A159) 


0 = -P _1 R/? 


(A160) 


Assuming 


0 = Be Xl (A161) 

then 

0 = ABe A ‘ (A162) 

Substituting (A161) and (A162) into (A160) yields a standard eigenvalue problem 


AB = -P -1 RB 


(A163) 


where A is an eigenvalue and B is an eigenvector. It should be noted that 


P" 1 = 


-M -1 C 

I 


M" 1 

0 


(A164) 


so that only M (whose dimension is half that of P) must be inverted. Many times 
in practice the interconnection forces defined by equations (A66) and (A67) contain no 
inertia terms. In these cases M becomes the identity and the inversion is trivial. 
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The eigenvalues and eigenvectors are, in general, complex quantities. The eigenvalues 
provide the natural frequencies and degree of damping for each mode in the system. 
Damping can be represented as a critical damping ratio for each mode using the relation 

Ci = - (A165) 

The eigenvectors (or modeshapes) give the relative shape of the deformation that oc- 
curs for a given mode. Since they are complex, displaying the vectors in a physically 
meaningful way is not entirely straightforward. For a given eigenvector B, it can be seen 
from equation (A156) that the lower half represents the displacement of the generalized 
coordinate, hence 


(A166) 

This part of the eigenvector (N), which is expressed in generalized coordinates, must be 
transformed back to physical coordinates using the transformations given in equations 
(A31 - (A34) and (A74) - (A78). Also, if condensation has been performed, the associated 
transformations must be reversed. This yields a complex eigenvector expressed in physical 
coordinates with the form 


Wj = { 


x 

y 

©z 

Z 

©y 

Xh 

yh 

©z h 

Zh 

©yh 


(A167) 


The eigenvectors and eigenvalues were used in equation (A101) to separate the variable (3 
into spatial and temporal factors. The transformations performed on (3 are now performed 


on B which, from equation (A161) yields 



w,<t) = w ie 


(A168) 


Kt 


This can be expanded in terms of real and imaginary parts as 
wj(t) = [Wr. + 

= [Wr. +;W I .]e AR - < [cos(A/ i t) + jsin(A/ i t)] 

= [Wr. cos(A /, 0 - W I; sinfA/.-Oje**.-* 

+ j[Wr. sin(A/,t) + W I; cos(A /t t)]e Afi .' (A169) 

Since w(t) must be real, only the real part is needed, so that the motion of the system 
due to a given mode is 


w;(t) = [Wr. cos(A u t) - W I; sin(A /i .t)]e A ' t .‘ (A170) 

The exponential term determines the rate of decay or growth of the motion for the 
mode. The modeshape characteristics can be displayed by plotting the motion through 
one period and neglecting the decay term. This will yield a three dimensional figure; 
however, the lateral motions (t/,-, 2 ,) at particular axial locations are usually of primary 
interest. For example, the rotor motion at the k th location for the i th mode will be 


Vk = Y Rkt cos(A lt t) - Y Ik , sin(A lt t) (A171) 

z k = Z Rki cos(A j.t) - Z Iki sin(A h t) (A172) 

To display these motions, it is necessary to let t vary such that Xj.t covers the range from 
0 to 27 t. The motions given by equations (A171) and (A172) will trace out an ellipse 
in the y — z plane. Housing motions and relative motions (rotor minus housing) can be 
displayed in an identical manner. It should be noted that the display is more than just 
the eigenvector. It is the transient motion of the vector (without the decay) through one 
period. 
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The forced response program obtains the particular solution to the equations of 
motion. The excitation must be harmonic and can be in one of three forms. The first 
and most often used form of excitation is due to mass unbalance of the rotor. This 
excitation is harmonic at a frequency equal to rotational speed. Since the force vector is 
rotating, the y and z axis components are ±90 degrees out of phase, depending on the 
direction of rotation. Since it is inertial, it is applied only to the rotor. The second form 
of excitation is static (zero frequency) and is due to circumferential pressure variations 
in the turbopump. It is applied to the rotor with an equal but opposite force applied 
to the housing. It is also distributed along the length of the rotor. The third form 
is due to rotating eccentricities in the interconnection elements. The magnitude of the 
force is equal to the product of the eccentricity and the element’s stiffness. For rotor 
eccentricities, the frequency is equal to rotation speed. For the rolling element bearings, 
the ball or roller train can become eccentric due to variations in element size and produce 
an eccentricity rotating at the speed of the rolling element separator (cage). The force is 
applied to the rotor with an equal but opposite force applied to the housing. The y and 
z axis components are ±90 degrees out of phase, depending on the direction of rotation. 
All three forms of excitation are discretely distributed along the length of the rotor. 

The equations of motion can be written for these types of excitation as 


M(fI)T7+C(f2)i7+K(fi)7?= Q(ut) (A173) 

where the dependency of M, C, and K on rotational speed (fi) has been emphasized and 

u = pQ. (A174) 

Depending on the form of excitation, p will be 1,0, or equal to the rolling element cage 
speed ratio. Since the equations are linear and a harmonic excitation is applied, the 
particular solution will be harmonic. The excitation can be written as 

Q(<) = K[Qe Jwt ] = »[(Qr + j'Qi)e ju,t ] (A175) 
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where 3?[ ] denotes the real part. The generalized force vector Q can be defined by 
comparison with the right hand side of equation (A84): 


Q<<) = ( ] ( A176 > 
l -*jFEh,M-*£F Eb ,(() J 

where F Ery , F Eri , F Ehy , and F Ehi , are the excitation forces for the rotor and housing 
along the y and z axes. For the unbalance excitation 

F E r v i{t) = (ma),fl 2 cos(Slt + <£,) 

(A177) 

F Er,i(t ) = (ma),fl 2 sin(f2t + <£,) 

(A178) 

Fkh„i(t) = F e h z i(t) = 0. 

(A179) 

For the static excitation 


F E r y = F Sy 

(A180) 

Ffir, = Fs, 

(A181) 

F E h y = F Ery 

(A182) 

f eu, = F Eri 

(A183) 

For the eccentricity excitation 


Fkr„«(0 = fc,€, COS (pilt + fa) 

(A184) 

F E r t i{t) = sin (pQt + <j>i) 

(A185) 

F Ehv i(t) = F E r t i(t) 

(A186) 
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F Bk , i(0 = F E r,i(t) (A187) 

From (A175) and (A176) 

F B (0 = »[(R + ;'I)e J ' pn< ] = R cos(pflt) - I sin(pQt) (A188) 

For unbalance excitation, from equations (A177) - (A178) 


R r i = ( ma),47 2 cos <f>i 

(A189) 

I r i = (ma),Sl 2 sin 4>i 

(A190) 

Rr t i = (ma),fl 2 sin </>,• 

(A191) 

I Tti = ~(ma)iQ.~ cos (f>i 

(A192) 

Rh y = Rh, = Ih v = h, = 0. 

(A193) 

For the static excitation, from equations (A180) - (A181) 


Rryi = F Sy , = |Fs, | cos tan- 1 

(A194) 

Rr,i = Fsz, = \Fs, 1 sin tan" 1 

(A195) 


Rh y i = Rr v i (A196) 

Rk.i = Rr,i (A197) 


and I T . and are of no consequence since in equation (A188) 

sin(pflt) = sin(Ot) = 0. (A198) 


For the eccentricity excitation, from equation (A184) - (A185) 
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Rr y i — COS (pi 

(A199) 

Ir v i = sin <pi 

(A200) 

Rr z i = sin <f>i 

(A201) 

Ir 2 i — ^ £ t COS <f>i 

(A202) 

Rk y i = Rr y i 

(A203) 

I hyi — 1 r y i 

(A204) 

Rfi 2 i — Rr z i 

(A205) 

= Ir z i 

(A206) 


From equations (A175), (A176), and (A188) 


(Qr + jQi) 


0 


4> T (R-,, + 3 Ir y ) 

* T ( R r. +yir.) 

+ A,) - *£( R h. + A.) 


(A207) 


where R and I are defined using either equations (A189) - (A193), (A194) - (A197), 
or (A199) - (A206) depending on the form of the excitation. With Q expressed as in 
equation (A175), the solution can be assumed to be of a similar form 


7?(t) = ft[Ne Ju;t ] 


(A208) 


so that 


and 


rj(t) = 


(A209) 
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f,(t) = ft[-w 2 Ne JU,t ] (A210) 

Substituting (A175) and (A208) - (A210) into (A173) and canceling the e ju)t term yields 


(K(Q) - u 2 M(fi) + ;wC(fi)]N = Q (A211) 

Equating real and imaginary parts of equation (A211) gives 

(K — w 2 M)Nr - wCNi = Qr (A212) 

wCNr + (K — w 2 M)Ni = Qi (A213) 

Equations (A212) and (A213) can be solved simultaneously for Nr and Ni. 

The response vector N defined by equation (A208) can be displayed in exactly the 
same way as the N defined as the displacement part of the eigenvector in equation (A166). 

The procedure is exactly the same except that the decay term that is ignored when 

displaying the eigenvector never appears for the response vector; therefore, the procedure 
will not be repeated here. 

The nonlinear transient simulation obtains the complete solution of equation (A84) 
with all of the nonlinearities and excitations that have been discussed. Rewriting equation 
(A84) with the introduction of rj as defined by equation (A155) gives 

Mrj + C17 + K77 = Q(t) (A214) 

This equation can be put into a form that is convenient for numerical solution. First, it 
can be rearranged as 


t 7 = M" 1 [Q(0-Kt 7 -Cfl 


(A215) 


If the interconnection forces defined by equation (A66) and (A67) contain no inertia 
terms, then M becomes the identity and inversion is not required. Equation (A215) can 
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be written with no derivatives on the right hand side by introducing a new variable and 
corresponding equation 


fj = v 


u = M- x [Q(t)-K»7-Ci/] 


(A216) 

(A217) 


This set of equations is now in the general form of 


x = f(x, t ) (A218) 

These equations will be integrated using the Adams-Moulton predictor-corrector method. 
The predictor equation for this method is 


x 


p 

n + l 


— x n + ^ [55f(x n ,« n ) 59f(x n 


The corrector equation is 


tn — i) "h 37f(x„_2, tn-?) (x n _3 , t n _ 3 )j 

(A219) 


X„+i = x n + ~ [9f(x£ +1> t n+ , )+ 19f(x n , t n ) - 5f(x n _i, t n _, ) + f(x n _ 2 , t n -a)] (A220) 

This method requires four starting points which can be obtained by using a simpler 
method. The method used here is the Euler or tangent line method 


x n +i = x „ + hf(x B ,t n ) (A221) 

The first point needed is the prescribed initial condition vector. Equation (A221) is used 
three times to obtain enough values to begin using equations (A219) and (A220). An 
alternate method is built into the simulation as an option. The alternate method was not 
used in this study and will not be described here. 

The numerical integration produces the generalized coordinate motion at the discrete 
time points 
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t n — n ■ h 


(A222) 


where n is an integer and varies from zero to the specified maximum. The generalized 
coordinate motion is transformed back into physical coordinates and these coordinates, 
along with interaction forces and housing accelerations can be displayed versus time. It is 
usually necessary to perform certain operations on the results after they have been gen- 
erated. This post-processing includes such things as filters and Fourier transformations. 

The M, C, and K matrices in equation (A217) are, in general, functions of engine 
power level. As discussed earlier, due to the relationship between power level and speed, 
they can be expressed as functions of speed. If speed is a function of time, these matrices 
(and M -1 if required) must be frequently re-evaluated. It is not necessary to re-evaluate 
them at each time step, however. The frequency for re-evaluation is specified in terms 
of a speed increment instead of time. This increment should be made small enough to 
keep the change in the interconnection coefficients small. This approach allows for more 
efficient numerical solution of the system equations. 
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APPENDIX B 

HPFTP MODEL NOMINAL DATA 


The nominal rotor-housing interconnection coefficient data axe shown in figures 113 
through 119. These data were provided by the SSME manufacturer with the exception of 
the pump interstage seal cross-coupled stiffness coefficients. The coefficients provided by 
the manufacturer were replaced with functions of the form of equation 125 where cr — .6. 
The resulting curve closely matched the original data. The damping coefficients for all 
bearings was a constant of 3.0 • The frequencies of the nine free- free rotor component 

modes are shown in table 2. The frequencies of the twenty free interface housing modes 
are shown in table 3. 


Table 2. Frequencies of free-free rotor modes. 

Mode 1 - 0.0 Hz. 

Mode 2 - 0.0 Hz. 

Mode 3 - 634.5 Hz. 

Mode 4 - 1350. Hz. 

Mode 5 - 1910. Hz. 

Mode 6 - 2591. Hz. 

Mode 7 - 3216. Hz. 

Mode 8 - 3935. Hz. 

Mode 9 - 3953. Hz. 
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Figure 115. HPFTP turbine interstage seal stiffness coefficient versus speed. 
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Figure 1 17. HPFTP turbine interstage seal damping coefficient versus speed. 




RPM 

Figure 118. HPFTP pump interstage seal cross-coupled stiffness coefficient versus speed. 





Table 3. Frequencies of free interface housing modes. 


Mode 1 - 50.29 Hz. 
Mode 2 - 114.0 Hz. 
Mode 3 - 363.9 Hz. 
Mode 4 - 417.9 Hz. 
Mode 5 - 712.9 Hz. 
Mode 6 - 836.3 Hz. 
Mode 7 - 920.1 Hz. 
Mode 8 - 995.8 Hz. 
Mode 9 - 1024. Hz. 
Mode 10 -1143. Hz. 
Mode 11 -1163. Hz. 
Mode 12 -1672. Hz. 
Mode 13 -1672. Hz. 
Mode 14 -1802. Hz. 
Mode 15 -1808. Hz. 
Mode 16 -2573. Hz. 
Mode 17 -2646. Hz. 
Mode 18 -2753. Hz. 
Mode 19 -3534. Hz. 
Mode 20 -3536. Hz. 
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